# This program is use for ENGAGE AFOLU paper
# 8, 03, 2020
# Contact person: Tomoko Hasegawa

library(gdxrrw)
library(ggplot2)
library(dplyr)
library(reshape2)
library(tidyr)
library(maps)
library(grid)
library(RColorBrewer)
library(scales)
library(cowplot)
library(car)
library(data.table)
#library(memisc)
library(tidyr)
library(plyr)

pastelpal1 <- brewer.pal(9, "Pastel1")
spectpal <- brewer.pal(11, "Spectral")
AR6pal <- c("C1: 1.5C with no or low OS"="#8fbc8f","C2: 1.5C with high OS"="#7fffd4","C3: lower 2C"="#6495ed","C4: higher 2C"="#f0e68c","C5: below 2.5C"="#ffa07a","C6: below 3.0C"="#ee82ee","C7: above 3.0C"="#a9a9a9")
variablemap <- read.table("../define/variablemap.txt",sep="\t",header=T) 
modelemap <- read.table("../define/modelmap.txt",sep="\t",header=T) 
dir.create("../output/")
dir.create("../output/fig/")
dir.create("../output/fig/CumCO2_single/")
dir.create("../output/fig/box/")
dir.create("../output/fig/box_multi/")
dir.create("../output/fig/area/")
dir.create("../output/fig/line_box/")
dir.create("../output/fig/line_box_FP/")
dir.create("../output/reg/")

MyThemeLine_grid <- theme_bw() +
  theme(
    panel.border=element_rect(fill=NA),
    panel.grid.minor = element_line(color = NA), 
    axis.line=element_line(colour="black"),
    panel.background=element_rect(fill = "white"),
    panel.grid.major=element_blank(),
    strip.background=element_rect(fill="white", colour="white"),
    strip.text.x = element_text(size=15, colour = "black", angle = 0,face="bold"),
    axis.text.x=element_text(angle=90, vjust=0.5, hjust=0.5),
    legend.text = element_text(size = 7),
    legend.title = element_text(size = 7)
  )

MyThemeLine_grid2 <- MyThemeLine_grid +
  theme(
    strip.text.x = element_text(size=10),
    legend.text = element_text(size = 10),
    legend.title = element_text(size = 10)
  )

MyThemeLine_grid2h <- MyThemeLine_grid2 +
  theme(
    axis.text.x=element_text(angle=0, vjust=0.5, hjust=0.5)
      )

colhead <- c("Model","Scenario","Region","Variable","Unit")
colhead2 <- c("ModelN","Scenario","Region","Variable2","Unit")
yearlist <- c("2005","2010","2015","2020","2025","2030","2035","2040","2045","2050","2055","2060","2065","2070","2075","2080","2085","2090","2095","2100")
yearxlist <- c("X2005","X2010","X2015","X2020","X2025","X2030","X2035","X2040","X2045","X2050","X2055","X2060","X2065","X2070","X2075","X2080","X2085","X2090","X2095","X2100")
#Read Data

#worlddata <- read.csv("../data/snapshot_world_iamc_engage_202003.csv", header=T) 
worlddata <- read.csv("../data/snapshot_world_iamc_engage_202009.csv", header=T) 

#worlddata <- worlddata[, c(-6,-7)] 
Modellist <- unique(worlddata$Model)
Scenariolist <- unique(worlddata$Scenario)
Variablelist <- as.vector(unique(worlddata$Variable))
worlddata1 <-worlddata %>% gather(key=Year,value=Value,-Model,-Scenario,-Region,-Variable,-Unit) %>% inner_join(variablemap,by="Variable")  %>%
  select(-Variable) %>% inner_join(modelemap,by="Model") %>% select(-Model)  %>% select(append(colhead2,c("Year","Value"))) 
row.has.na.worlddata1 <- apply(worlddata1, 1, function(x){any(is.na(x))})
worlddata1 <- worlddata1[!row.has.na.worlddata1,] 

#-----Regional data
scenariolist_engage <- read.table("../individual/ENGAGE/scenariolist_engage_plot2.txt",sep="\t",header=T) 
regionaldata <- read.csv("../data/snapshot_Compare_R5.csv", header=T) 
#regionaldata <- read.csv("../data/snapshot_Compare_R10.csv", header=T) 
varreg <- read.table("../individual/ENGAGE/VAR_region.txt", header=T, sep="\t") 
regionmap <- read.table("../individual/ENGAGE/Region_R5.map", header=T, sep="\t") 
#regionmap <- read.table("../individual/ENGAGE/Region_R10.map", header=T, sep="\t") 

regionaldata1 <-regionaldata %>% gather(key=Year,value=Value,-Model,-Scenario,-Region,-Variable,-Unit) %>% inner_join(variablemap,by="Variable")  %>%
  select(-Variable) %>% inner_join(modelemap,by="Model") %>% select(-Model) %>% 
  inner_join(varreg,by=c("Variable2")) %>% inner_join(scenariolist_engage,by=c("Scenario"))  %>% select(append(colhead2,c("Year","Value"))) 
row.has.na.regionaldata1 <- apply(regionaldata1, 1, function(x){any(is.na(x))})
regionaldata1 <- regionaldata1[!row.has.na.regionaldata1,]
regionaldata1.0 <- regionaldata1 %>% inner_join(regionmap,by=c("Region")) %>% select(-Region)
regionaldata1.0 <- regionaldata1.0[c(1,2,7,3,4,5,6)]
names(regionaldata1.0) <- c("ModelN","Scenario","Region","Variable2","Unit","Year","Value")	#Rename the column
worlddata1.0 <- rbind(worlddata1,regionaldata1.0)

#-----Data preparation
variableplot1_load <- read.table("../define/variablelist1.txt", header=F, sep=",") 
variableplot2_load <- read.table("../define/variablelist2.txt", header=F, sep=",") 
variableplot_load <- rbind(variableplot1_load,variableplot2_load)
variableplot <- as.vector(variableplot_load$V1)
variableplot_load <- read.table("../define/variablelist_engage2.txt", header=F, sep=",") 
variableplot <- as.vector(variableplot_load$V1)

linepalette <- c("#4DAF4A","#377EB8","#E41A1C","#FF7F00","#984EA3","#FFFF33","#A65628","#F781BF","#8DD3C7","#FB8072","#80B1D3","#FDB462","#B3DE69","#FCCDE5","#D9D9D9","#BC80BD","#CCEBC5","#FFED6F")

#CO2 emissions and cumulative emissions
worldCO2 <- worlddata1.0 %>% filter(Variable2=="Emi_CO2")
row.has.na.co2 <- apply(worldCO2, 1, function(x){any(is.na(x))})
worldCO21 <- worldCO2[!row.has.na.co2,] 
CarPri <- worlddata1 %>% filter(Variable2=="Price_Carbon")
row.has.na.CarPri <- apply(CarPri, 1, function(x){any(is.na(x))})
CarPri1 <- CarPri[!row.has.na.CarPri,] 
symDim <- 7
attr(worldCO21, "symName") <- "IPCCSR15WorldCO2Emi"
attr(CarPri1, "symName") <- "IPCCSR15WorldCarPrice"
lst2 <- wgdx.reshape(worldCO21, symDim)
wgdx.lst(gdxName = "../output/WorldCO2all.gdx",lst2)
lst3 <- wgdx.reshape(CarPri1, symDim)
wgdx.lst(gdxName = "../output/CarbonPrice.gdx",lst3)

attr(worlddata1.0, "symName") <- "IPCCSR15All"
lst4 <- wgdx.reshape(worlddata1.0, symDim)
wgdx.lst(gdxName = "../output/IPCCSR15All.gdx",lst4)

system(paste("gams CO2.gms",sep=" "))
CumEmi <- rgdx.param('../output/cumulativeEmissions.gdx','CumCO2Emis')
TCRECumEmi <- rgdx.param('../output/cumulativeEmissions.gdx','TCRECumCO2Emis')
TCRECumEmi$TCRECumCO2Emis[TCRECumEmi$TCRECumCO2Emis>5000] <- 4900
PeakWarming <- rgdx.param('../output/IPCCSR15allPlus.gdx','PeakWarmingF')
names(PeakWarming) <- c("ModelN","Scenario","WPeak","PeakWarmingF")
Warming2100 <- rgdx.param('../output/IPCCSR15allPlus.gdx','Warming2100F')
names(Warming2100) <- c("ModelN","Scenario","W2100","Warming2100F")
WarmingF <- inner_join(PeakWarming,Warming2100) %>% select(-PeakWarmingF,-Warming2100F)

worlddata1Plus <- rgdx.param('../output/IPCCSR15allPlus.gdx','IPCCSR15AllPlus')
names(worlddata1Plus) <- c("ModelN","Scenario","Region","Variable2","Unit","Year","Value")
worlddata1.0 <- worlddata1Plus
Worlddata1.1_CumCO2 <- left_join(worlddata1.0,select(TCRECumEmi,-Unit),by=c("ModelN","Scenario","Region"))
Unitlist <- unique(select(worlddata1.0,c("Variable2","Unit")))
Unitlist2 <- left_join(dplyr::rename(variableplot_load,"Variable2"=V1),Unitlist)

NetZeroYear <- rgdx.param('../output/cumulativeEmissions.gdx','NetZeroTimingCO2Emis')

NetZeroYear <- NetZeroYear %>% spread(key=Year, value=NetZeroTimingCO2Emis, fill=NA)
#names(NetZeroYear) <- c("ModelN","Scenario","Region","Variable3","Unit","1990","2005","2025","2030","2035","2040","2045","2050","2055","2060","2065","2070","2075","2080","2085","2090","2095","2100","2110")	#Rename the column
names(NetZeroYear) <- c("ModelN","Scenario","Region","Variable3","Unit","2020","2025","2030","2035","2040","2045","2050","2055","2060","2065","2070","2075","2080","2085","2090","2095","2100","2110")	#Rename the column
# remove 2110
NetZeroYear <- NetZeroYear[c(-23)]
NetZeroYear <- NetZeroYear %>% gather(key=NetZeroYear,value=Value,-ModelN,-Scenario,-Region,-Variable3,-Unit)
NetZeroYear <- na.omit(NetZeroYear)
NetZeroYear$NetZeroYear <- as.numeric(NetZeroYear$NetZeroYear)

Worlddata1.1_CumCO2_NetZero <- left_join(Worlddata1.1_CumCO2,select(NetZeroYear,-Unit,-Variable3,-Value),by=c("ModelN","Scenario","Region"))

AFOLUNetZeroYear <- rgdx.param('../output/cumulativeEmissions.gdx','NetZeroTimingAFOLUEmis')
AFOLUNetZeroYear <- AFOLUNetZeroYear %>% spread(key=Year, value=NetZeroTimingAFOLUEmis, fill=NA)
#names(AFOLUNetZeroYear) <- c("ModelN","Scenario","Region","Variable3","Unit","2010","2025","2030","2035","2040","2045","2050","2055","2060","2065","2070","2075","2080","2085","2090","2095","2100","2110")	#Rename the column
names(AFOLUNetZeroYear) <- c("ModelN","Scenario","Region","Variable3","Unit","2005","2025","2030","2035","2040","2045","2050","2055","2060","2065","2070","2075","2080","2085","2090","2095","2100")	#Rename the column
AFOLUNetZeroYear <- AFOLUNetZeroYear %>% gather(key=AFOLUNetZeroYear,value=Value,-ModelN,-Scenario,-Region,-Variable3,-Unit)
AFOLUNetZeroYear <- na.omit(AFOLUNetZeroYear)
AFOLUNetZeroYear$AFOLUNetZeroYear <- as.numeric(AFOLUNetZeroYear$AFOLUNetZeroYear)

#-----Engage figure
dir.create("../output/fig/ENGAGE/")
dir.create("../output/fig/ENGAGE/box/")
dir.create("../output/fig/ENGAGE/line/")
dir.create("../output/fig/ENGAGE/schatter/")
dir.create("../output/fig/ENGAGE/txt")
dir.create("../output/fig/ENGAGE/line_box/")

#Extract Engage data
scecbfpmap <- read.table("../define/scecbfpmap.txt",sep="\t",header=T) 
Engagepal <- c("Full"="red","Peak"="blue")
shapepallet <- c('AIM_CGE 2_2'=16, 'COFFEE 1_1'=17,'IMAGE 3_0'=11, 'MESSAGEix-GLOBIOM_1.0'=15, 'POLES_ENGAGE'=3, 'REMIND-MAgPIE 2_0-4_1'=7, 'WITCH 5_0'=8, 'GEM-E3_092019'=2)
fillpalland <- c('Pasture'='lightblue','Forest'='lightgreen','Cropland for Energy'='mediumpurple','Cropland for NonEnergy'='salmon')
varmodel_unused <- read.table("../individual/ENGAGE/VAR_Model_unused.map",header=T, sep = "\t")
model_unused <- read.table("../individual/ENGAGE/Model_unused.map",header=T, sep = "\t")
scenariomodel_unused <- read.table("../individual/ENGAGE/Scenario_Model_unused.map",header=T, sep = "\t")


worlddata3.0 <- worlddata1.0 %>% inner_join(scecbfpmap,by="Scenario") %>% inner_join(scenariolist_engage,by=c("Scenario")) %>%
  anti_join(varmodel_unused,by=c("Variable2","ModelN")) 
worlddata3.0$Year <- as.numeric(levels(worlddata3.0$Year))[worlddata3.0$Year]
worlddata3.0$CarbonBudget <- as.numeric(levels(worlddata3.0$CarbonBudget))[worlddata3.0$CarbonBudget]
#worlddata3.0$CarbonBudget <- as.numeric(worlddata3.0$CarbonBudget)
worlddata3.0 <- worlddata1.0 %>% inner_join(scecbfpmap,by="Scenario") %>% inner_join(scenariolist_engage,by=c("Scenario")) %>%
  anti_join(varmodel_unused,by=c("Variable2","ModelN")) 
worlddata3.0$Year <- as.numeric(levels(worlddata3.0$Year))[worlddata3.0$Year]
worlddata3.0$CarbonBudget <- as.numeric(levels(worlddata3.0$CarbonBudget))[worlddata3.0$CarbonBudget]
#Scatter figure
regresults <- as.list(variableplot)
names(regresults) <- names(variableplot)
regresults2 <- as.list(variableplot)
names(regresults2) <- names(variableplot)
figurelist <- as.list(variableplot)
names(figurelist) <- names(variableplot)
#y51 <- c(2030,2050)
y51 <- c(2050,2090)

yearmap <- read.table("../define/yearmap.txt",sep="\t",header=T) 

# calculate mean values
worlddata3.0.1 <- worlddata3.0 %>% anti_join(scenariomodel_unused,by=c("Scenario","ModelN"))
worlddata3.0.1 <- filter(worlddata3.0.1, Region %in% c("World"), Year %in% c(2040, 2050, 2060, 2080, 2090, 2100) & CarbonBudget<=2000)
worlddata3.0.1 <-worlddata3.0.1 %>% inner_join(yearmap,by="Year")
worlddata3.0.2 <- ddply(worlddata3.0.1,c("Variable2","Term","Full_Peak"),transform,v50=quantile(Value,0.50))
worlddata3.0.2mean <- worlddata3.0.2[worlddata3.0.2$ModelN=="MESSAGEix-GLOBIOM_1.0", c(-1)]

worlddata3.0.3 <- ddply(worlddata3.0.2,c("Variable2","Year","Full_Peak"),transform,v50y=quantile(Value,0.50))

for (i in 1:length(variableplot)){
  worlddata3.1 <- filter(worlddata3.0, Region %in% c("World"),Variable2==variableplot[i])
  worlddata3.1 <- worlddata3.1 %>% anti_join(scenariomodel_unused,by=c("Scenario","ModelN"))
  varname <- variableplot[i]
  if(nrow(worlddata3.1)>=2){
    worlddata3.2 <- filter(worlddata3.1,Year %in% c(2050) & CarbonBudget<=2000)
#    worlddata3.2 <- filter(worlddata3.1,Year %in% c(2040, 2045, 2050, 2055, 2060) & CarbonBudget<=2000)
#    worlddata3.2 <- filter(worlddata3.1,Year %in% c(2040, 2050, 2060, 2080, 2090, 2100) & CarbonBudget<=2000)
    worlddata3.2 <-worlddata3.2 %>% inner_join(yearmap,by="Year")
    worlddata3.2.0 <- worlddata3.2 %>% filter(CarbonBudget!=300 & CarbonBudget!=400 & CarbonBudget!=500 & CarbonBudget!=1600 & CarbonBudget!=2000)
    worlddata3.2.1 <- ddply(worlddata3.2.0,c("CarbonBudget","Term","Full_Peak"),transform,min=quantile(Value,0.33))
    worlddata3.2.2 <- ddply(worlddata3.2.1,c("CarbonBudget","Term","Full_Peak"),transform,max=quantile(Value,0.66))
    worlddata3.2.3 <- ddply(worlddata3.2.2,c("CarbonBudget","Term","Full_Peak"),transform,v50=quantile(Value,0.50))
    
    maxv <- max(worlddata3.2$Value)
    minv <- min(worlddata3.2$Value)
    
    gschatter <- ggplot(worlddata3.2,aes(x=CarbonBudget, y = Value, colour=Full_Peak)) +  
#      geom_boxplot(aes(x=CarbonBudget, y = Value, fill=Full_Peak), size=0.5,alpha=0.7,outlier.color = NA) + 
      geom_point(aes(shape=ModelN), size=2, alpha=0.7,outlier.color = NA) + 
      geom_smooth(method='lm', se = FALSE, size=0.5) +
      geom_ribbon(data=worlddata3.2.3,aes(x=CarbonBudget,ymin=min, ymax=max, fill=Full_Peak), alpha=0.1) +
    #  geom_line(data=worlddata3.2.3,aes(x=CarbonBudget,y=v50), size = 0.5, colour="black") +
      annotate("rect",xmin=350,xmax=770,ymin=minv,ymax=maxv,alpha=0.2,fill=spectpal[10]) +
      annotate("rect",xmin=1250,xmax=1690,ymin=minv,ymax=maxv,alpha=0.2,fill=spectpal[8]) +
      ylab(Unitlist2$Unit[i]) + 
      #ylab(labelname(variableplot[i])) +
      xlab("Carbon budget [GtCO2]") +labs(fill="")+ MyThemeLine_grid + 
      scale_colour_manual(name="Budget cap",values=Engagepal) +
      scale_fill_manual(name="Budget cap",values=Engagepal) +
      scale_shape_manual(values=shapepallet) +
      theme(text=element_text(size=12),  
            axis.text.x=element_text(angle=0, vjust=0.9, hjust=0.5, size = 12),
            #legend.position="none"
            ) +
      ggtitle(label=variableplot[i]) +
      guides(color=guide_legend(ncol=1))
#    figurelist[[variableplot[i]]] <- gschatter +facet_wrap(Term~.,ncol=4)
    outname <- paste0("../output/fig/ENGAGE/schatter/",variableplot[i],".png")
#    outname <- paste0("../output/fig/ENGAGE/schatter_poles/",variableplot[i],".png")
    ggsave(gschatter, file=outname, dpi = 600, width=6, height=4,limitsize=FALSE)

#    gschatter2 <- gschatter +facet_wrap(Term~.,ncol=4)
#   plot(gschatter)
  }
}

for (i in 1:length(variableplot)){
  worlddata3.1 <- filter(worlddata3.0, Region %in% c("World"),Variable2==variableplot[i])
  worlddata3.1 <- worlddata3.1 %>% anti_join(scenariomodel_unused,by=c("Scenario","ModelN"))
  varname <- variableplot[i]
  if(nrow(worlddata3.1)>=2){
    if(length(as.vector(unique(worlddata3.1$ModelN)))>=2){
      for(k in 1:length(y51)){

        #Regression
        worlddata3.3 <- worlddata3.1 %>% filter(Year %in% c(y51[k]-10,y51[k],y51[k]+10)  & CarbonBudget<=2000)
        regression <- lm(data=worlddata3.3,Value ~ CarbonBudget+Full_Peak+ModelN)
        regression2 <- lm(data=worlddata3.3,Value ~ Full_Peak*CarbonBudget+ModelN)
        regresults[[variableplot[i]]] <- summary(regression)
        regresults2[[variableplot[i]]] <- summary(regression2)
        
        if (i+k==2){
          reg_summ <- cbind(as.character(variableplot[i]),y51[k],tibble::rownames_to_column(as.data.frame(regresults[[variableplot[i]]]$coefficients), "factors"))
          reg_summ2 <- cbind(as.character(variableplot[i]),y51[k],tibble::rownames_to_column(as.data.frame(regresults2[[variableplot[i]]]$coefficients), "factors"))
          num_obs <- cbind(as.character(variableplot[i]), nrow(worlddata3.3))
        }else{
          reg_summ <- rbind(reg_summ, cbind(as.character(variableplot[i]),y51[k],tibble::rownames_to_column(as.data.frame(regresults[[variableplot[i]]]$coefficients), "factors")) )
          reg_summ2 <- rbind(reg_summ2, cbind(as.character(variableplot[i]),y51[k],tibble::rownames_to_column(as.data.frame(regresults2[[variableplot[i]]]$coefficients), "factors")) )
          num_obs <- rbind(num_obs,cbind(as.character(variableplot[i]), nrow(worlddata3.3)))
        }

        # Residual analysis
        outname <- paste0("../output/fig/ENGAGE/regression/analysis/",variableplot[i],".png")
        par(mfrow=c(2,2))
        plot(regression)
        dev.print(png, file = outname, width = 1000, height = 800)
        
        pairsdf <- worlddata3.3[c(1,7,8,9)]
        pairs(pairsdf)
        outname <- paste0("../output/fig/ENGAGE/regression/pairs/",variableplot[i],".png")
        dev.print(png, file = outname, width = 800, height = 600)
        
      }
    }
  }
}

write.table(reg_summ,file="../output/fig/ENGAGE/regression/regression.txt", append = FALSE, row.names=TRUE, quote = FALSE, sep = ",")
write.table(reg_summ2,file="../output/fig/ENGAGE/regression/regression2.txt", append = FALSE, row.names=TRUE, quote = FALSE, sep = ",")
reg_summ_df <- as.data.frame(reg_summ)
reg_summ2_df <- as.data.frame(reg_summ2)
num_obs_df <- as.data.frame(num_obs)

names(reg_summ_df) <- c("factor","Year","Indicators","coefficient","Std. Error","t-value","Pr>t")
names(reg_summ2_df) <- c("factor","Year","Indicators","coefficient","Std. Error","t-value","Pr>t")
names(num_obs_df) <- c("factor","Num_obs")

reg_summ_df$coefficient <- as.numeric(reg_summ_df$coefficient)
reg_summ_df$`Std. Error` <- as.numeric(reg_summ_df$`Std. Error`)
reg_summ_df$`t-value` <- as.numeric(reg_summ_df$`t-value`)
reg_summ_df$`Pr>t` <- as.numeric(reg_summ_df$`Pr>t`)
reg_summ2_df$coefficient <- as.numeric(reg_summ2_df$coefficient)
reg_summ2_df$`Std. Error` <- as.numeric(reg_summ2_df$`Std. Error`)
reg_summ2_df$`t-value` <- as.numeric(reg_summ2_df$`t-value`)
reg_summ2_df$`Pr>t` <- as.numeric(reg_summ2_df$`Pr>t`)
#num_obs_df$`Num_obs` <- as.numeric(num_obs_df$`Num_obs`)

#reg_summ_df <- tibble::rownames_to_column(reg_summ_df, "factors")
#reg_summ2_df <- tibble::rownames_to_column(reg_summ2_df, "factors")
#tmp <- tibble::rownames_to_column(as.data.frame(regresults[[variableplot[i]]]$coefficients), "factors")

symDim <- 4
symDim1 <- 2
attr(reg_summ_df, "symName") <- "regression1"
attr(reg_summ2_df, "symName") <- "regression2"
attr(num_obs_df, "symName") <- "num_obs"
lst2 <- wgdx.reshape(reg_summ_df, symDim)
lst3 <- wgdx.reshape(reg_summ2_df, symDim)
lst <- wgdx.reshape(num_obs_df, symDim1)
wgdx.lst(gdxName = "../output/fig/ENGAGE/regression/reg1.gdx",lst2)
wgdx.lst(gdxName = "../output/fig/ENGAGE/regression/reg2.gdx",lst3)
wgdx.lst(gdxName = "../output/fig/ENGAGE/regression/num_obs.gdx",lst)

#------Residual analysis (uncompleted!!!)
system(paste("gams regression_analysis.gms",sep=" "))

#Residual plots
moddummyave <- rgdx.param
fpmap <- read.table("../define/fpmap.txt",sep="\t",header=T) 
lnormframe <- 0
for (i in 1:length(variableplot)){
  worlddata3.1 <- filter(worlddata3.0, Region %in% c("World"),Variable2==variableplot[i])
  varname <- variableplot[i]
  if(nrow(worlddata3.1)>=2){
    if(length(as.vector(unique(worlddata3.1$ModelN)))>=2){
      for(k in 1:length(y51)){
        worlddata3.3 <- worlddata3.1 %>% filter(Year %in% c(y51[k]-10,y51[k],y51[k]+10)  & CarbonBudget<=2000)
        worlddata3.3 <- worlddata3.3 %>% inner_join(fpmap,by="Full_Peak")

        #regression <- lm(data=worlddata3.3,Value ~ CarbonBudget+Full_Peak+ModelN)
        
        CBela <- reg_summ_df[reg_summ_df$factor==varname & reg_summ_df$Indicators=="CarbonBudget" & reg_summ_df$Year==y51[k],c(4)]
        FPela <- reg_summ_df[reg_summ_df$factor==varname & reg_summ_df$Indicators=="Full_PeakPeak" & reg_summ_df$Year==y51[k],c(4)]
        intercept <- reg_summ_df[reg_summ_df$factor==varname & reg_summ_df$Indicators=="(Intercept)" & reg_summ_df$Year==y51[k],c(4)] + moddummyave[moddummyave$var==varname & moddummyave$year==y51[k],c(3)] 


        CBval  <- worlddata3.3[worlddata3.3$Variable2==varname & worlddata3.3$Region=="World", c(8)]
        FPval  <- worlddata3.3[worlddata3.3$Variable2==varname & worlddata3.3$Region=="World", c(10)]
        obsval <- worlddata3.3[worlddata3.3$Variable2==varname & worlddata3.3$Region=="World", c(7)]
        
        estval <- CBval*CBela + FPval*FPela + intercept
        
        residual <- obsval - estval 

        lnormframe0 <- data.frame(Variable=variableplot[i],Year=y51[k],residual)
        
        if (i==1 & k==1){
          lnormframe <- lnormframe0
        }else{
          lnormframe <- rbind(lnormframe,lnormframe0)
        }
        
        
      }
    }
  }
}


#-----------------
scecbfpmap <- read.table("../define/scecbfpmap.txt",sep="\t",header=T) 
Engagepal <- c("Full"="red","Peak"="blue")
shapepallet <- c('AIM_CGE 2_2'=16, 'COFFEE 1_1'=17,'IMAGE 3_0'=11, 'MESSAGEix-GLOBIOM_1.0'=15, 'POLES_ENGAGE'=3, 'REMIND-MAgPIE 2_0-4_1'=7, 'WITCH 5_0'=8, 'GEM-E3_092019'=2)
fillpalland <- c('Pasture'='lightblue','Forest'='lightgreen','Cropland for Energy'='mediumpurple','Cropland for NonEnergy'='salmon')
varmodel_unused <- read.table("../individual/ENGAGE/VAR_Model_unused.map",header=T, sep = "\t")
model_unused <- read.table("../individual/ENGAGE/Model_unused.map",header=T, sep = "\t")
scenariomodel_unused <- read.table("../individual/ENGAGE/Scenario_Model_unused.map",header=T, sep = "\t")

worlddata3.0 <- worlddata1.0 %>% inner_join(scecbfpmap,by="Scenario") %>% inner_join(scenariolist_engage,by=c("Scenario")) %>%
  anti_join(varmodel_unused,by=c("Variable2","ModelN")) 

#Making multi figures
filename_land <- paste('../output/fig/ENGAGE/schatter/land.png', sep = '')
filename_emit <- paste('../output/fig/ENGAGE/schatter/emissions.png', sep = '')
filename_prodcons <- paste('../output/fig/ENGAGE/schatter/prodcons.png', sep = '')

#7 in 1 figure
pp1 <- plot_grid(figurelist[["Land_Cover_Cropland_NonEnergy_Crops_Change_2010"]],figurelist[["Land_Cover_Cropland_Energy_Crops_Change_2010"]],figurelist[["Land_Cover_Pasture_Change_2010"]],figurelist[["Land_Cover_Forest_Afforestation_and_Reforestation_Change_2010"]],figurelist[["Land_Cover_Forest_Natural_Change_2010"]],figurelist[["Land_Cover_Forest_Managed_Change_2010"]],figurelist[["Land_Cover_Forest_Change_2010"]],ncol=2)
ggsave(pp1, file=filename_land, dpi = 250, width=14, height=14,limitsize=FALSE)

pp2 <- plot_grid(figurelist[["Food_Demand"]],figurelist[["Food_Demand_Crops"]],figurelist[["Food_Demand_Livestock"]],figurelist[["Agricultural_Pro_Non-Energy_Crops"]],figurelist[["Agricultural_Pro_Non-Energy_Livestock"]],figurelist[["Agricultural_Pro_Energy_Crops"]],figurelist[["Primary_Energy_Biomass_Modern"]],ncol=2)
ggsave(pp2, file=filename_prodcons, dpi = 250, width=14, height=14,limitsize=FALSE)

#4 in 1 figure
pp3 <- plot_grid(figurelist[["Emi_CH4_N2O_AFOLU"]],figurelist[["Emi_CO2_AFOLU"]],figurelist[["Carbon_Sequestration_CCS_Biomass"]],figurelist[["Carbon_Sequestration_Land_Use"]],ncol=1)
ggsave(pp3, file=filename_emit, dpi = 250, width=11, height=14,limitsize=FALSE)


#Time series figure
scenariolist_engage3 <- read.table("../individual/ENGAGE/scenariolist_engage_plot3.txt",sep="\t",header=T) 
worlddata3.3 <- worlddata1.0
worlddata3.3$Year <- as.numeric(levels(worlddata3.3$Year))[worlddata3.3$Year]
######
worlddata3.3 <- worlddata3.0
######
worlddata3.4 <- left_join(worlddata3.3,WarmingF)
interact <- interaction(worlddata3.4$Scenario,worlddata3.4$ModelN, sep=" : ")
interact2 <- interaction(worlddata3.4$W2100,worlddata3.4$Full_Peak, sep=" : ")
worlddata3.5 <-cbind(worlddata3.4,interact,interact2)
#worlddata3.5 <- left_join(worlddata3.4,WarmingF)
worlddata3.5 <- na.omit(worlddata3.5)
worlddata3.5$W2100 <- factor(worlddata3.5$W2100, levels=c('Baseline','Over2.5Deg','2.5Deg','2Deg','1.5Deg'))
# Time series figure FvsP
worlddata3.5.1 <- filter(worlddata3.5,W2100 %in% c("2Deg","1.5Deg"))
worlddata3.5.1$interact2 <- factor(worlddata3.5.1$interact2, levels=c('2Deg : Full','1.5Deg : Full','2Deg : Peak','1.5Deg : Peak'))

linepalette <- c("#4DAF4A","#377EB8","#E41A1C","#FF7F00","#984EA3","#A65628","#F781BF","#FFFF33","#8DD3C7","#FB8072","#80B1D3","#FDB462","#B3DE69","#FCCDE5","#D9D9D9","#BC80BD","#CCEBC5","#FFED6F")
scenariopalette <- linepalette
scenariopalette_low <- c('Over2.5Deg'='grey60', '2.5Deg'='grey80', '2Deg'='green', '1.5Deg'='blue')
scenariopalette_low2 <- c('1.5Deg : Full'='red', '2Deg : Full'='orange', '1.5Deg : Peak'='blue', '2Deg : Peak'='green')
scenariopalette_low3 <- c("EN_NPi2020_600f"='red',"EN_NPi2020_600"='blue')


for (i in 1:length(variableplot)){
  worlddata3.6 <- filter(worlddata3.5,Variable2==variableplot[i]) 
  worlddata3.6.0 <- filter(worlddata3.6,Year %in% c(2000, 2010, 2020, 2030, 2040, 2050, 2060, 2070, 2080, 2090, 2100))
  worlddata3.6.1 <- ddply(worlddata3.6.0,c("W2100","Year"),transform,min=quantile(Value,0.33))
  worlddata3.6.2 <- ddply(worlddata3.6.1,c("W2100","Year"),transform,max=quantile(Value,0.66))
  worlddata3.6.3 <- ddply(worlddata3.6.2,c("W2100","Year"),transform,v50=quantile(Value,0.50))

  if(nrow(worlddata3.6)>=2){
    minv <- min(worlddata3.6[7])
    maxv <- max(worlddata3.6[7])
    
    gline <- ggplot() +  
      geom_line(data=worlddata3.6.3, aes(x=Year, y = v50, group=interact, colour=W2100), size=1,alpha=0.7) + 
#      geom_point(data=worlddata3.6, aes(x=Year, y= Value, colour=W2100, shape=ModelN),size=0.8) +
      geom_ribbon(data=worlddata3.6.3,aes(x=Year,ymin=min, ymax=max, fill=W2100), alpha=0.1) +
      ylab(Unitlist2$Unit[i]) + xlab("Years") +labs(fill="")+ MyThemeLine_grid + 
      ylim(minv-abs(minv)*0.1, maxv+abs(maxv)*0) +
      scale_shape_manual(values=shapepallet) +
      scale_fill_manual(name="Scenario",values=scenariopalette_low) +
      scale_colour_manual(name="Scenario",values=scenariopalette_low) + 
            theme(legend.position="none", text=element_text(size=12),  
            axis.text.x=element_text(angle=0, vjust=0.9, hjust=1, size = 12)) +
      ggtitle(label=variableplot[i]) +
      guides(color=guide_legend(ncol=1))
#    plot(gline)

    worlddata3.7 <- filter(worlddata3.6,Year %in% c(2050,2100))
    gbox <- ggplot() +  
      geom_boxplot(data=worlddata3.7,aes(x=W2100, y = Value , fill=W2100), size=0.5,alpha=0.7,range=1) +
      ylab(Unitlist2$Unit[i]) + xlab("Scenarios") +labs(fill="")+ MyThemeLine_grid +
      scale_fill_manual(name="Scenario",values=scenariopalette_low) +
      ggtitle(label=variableplot[i]) +
      guides(color=guide_legend(ncol=1)) +
      ylim(minv-abs(minv)*0.1, maxv+abs(maxv)*0) +
      ylab("") +
      xlab("") +
      theme(legend.position="none",text=element_text(size=12),  
            axis.text.x=element_text(angle=45, vjust=0.9, hjust=1, size = 13)) +
      theme(
        legend.position = "none",
        axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.line = element_line(colour = "white"),
        plot.title = element_text(colour = "white"),
        panel.border = element_blank(),
        axis.ticks=element_blank()
      )
    
    gbox2 <- gbox +facet_wrap(Year~.,ncol=4,strip.position = "bottom") + 
             theme(strip.text.x = element_text(size =10))

#    plot(gbox2)
    outname3 <- paste0("../output/fig/ENGAGE/line_box/",variableplot[i],".png")  
    multfigure2(gline,gbox2,outname3)
    
  }
}


multfigure2 <- function(p1,p2,x){
  grid.newpage()
  l <- grid.layout(nrow=1, ncol=2, heights=unit(c(1, 1.25),c("null","null")), widths=unit(c(1.9, 0.50),c("null","null")))
  png(filename=x, width=1000, height=800,res=200)
  v <- viewport(layout = l)
  pushViewport(v)
  print(p1, vp = viewport(layout.pos.row=1,layout.pos.col=1))
  print(p2, vp = viewport(layout.pos.row=1,layout.pos.col=2))
  popViewport()
  dev.off()  
  return()
}


# FvsP
for (i in 1:length(variableplot)){
  worlddata3.6 <- filter(worlddata3.5.1,Variable2==variableplot[i])
  worlddata3.6.0 <- filter(worlddata3.6,Year %in% c(2000, 2010, 2020, 2030, 2040, 2050, 2060, 2070, 2080, 2090, 2100))
  worlddata3.6.1 <- ddply(worlddata3.6.0,c("W2100","Year","Full_Peak"),transform,min=quantile(Value,0.33))
  worlddata3.6.2 <- ddply(worlddata3.6.1,c("W2100","Year","Full_Peak"),transform,max=quantile(Value,0.66))
  worlddata3.6.3 <- ddply(worlddata3.6.2,c("W2100","Year","Full_Peak"),transform,v50=quantile(Value,0.50))

  
  if(nrow(worlddata3.6)>=2){
    minv <- min(worlddata3.6[7])
    maxv <- max(worlddata3.6[7])
    
    gline <- ggplot(data=worlddata3.6.3) +  
      geom_ribbon(aes(x=Year,ymin=min, ymax=max, fill=interact2), alpha=0.2) +
      #      geom_line(data=worlddata3.6, aes(x=Year, y = Value, group=interact, colour=interact2), size=0.1,alpha=0.7) + 
      geom_line(aes(x=Year, y = v50, group=interact, colour=interact2), size=1,alpha=0.7) + 
      #      geom_line(data=worlddata3.6.3, aes(x=Year, y = v50, group=interact, colour=W2100), size=0.1,alpha=0.7) + 
      geom_point(data=worlddata3.6, aes(x=Year, y= Value, colour=interact2, shape=ModelN),size=2) +
      ylab(Unitlist2$Unit[i]) + xlab("Years") +labs(fill="")+ MyThemeLine_grid2 + 
      ylim(minv-abs(minv)*0.1, maxv+abs(maxv)*0) +
      scale_shape_manual(values=shapepallet) +
      scale_fill_manual(name="Scenario",values=scenariopalette_low2) +
      scale_colour_manual(name="Scenario",values=scenariopalette_low2) + 
      theme(
           legend.position="none", 
           text=element_text(size=14),  
           plot.title = element_text(size=12),
           axis.text.x=element_text(angle=0, vjust=0.9, hjust=0.5, size = 12)) +
      ggtitle(label=variableplot[i]) +
      guides(color=guide_legend(ncol=1))
   plot(gline)
#    outname <- paste0("../output/fig/ENGAGE/line_box_FP/legend.png")  
#    ggsave(gline, file=outname, dpi = 600, width=9, height=5,limitsize=FALSE)

    worlddata3.7 <- filter(worlddata3.6,Year %in% c(2050,2100))
    gbox <- ggplot() +  
      #      geom_boxplot(data=worlddata3.7,aes(x=W2100, y = Value , fill=W2100), size=0.5,alpha=0.7,range=1) +
      geom_boxplot(data=worlddata3.7,aes(x=interact2, y = Value , fill=interact2), size=0.5,alpha=0.7) +
      ylab(Unitlist2$Unit[i]) + xlab("Scenarios") +labs(fill="")+ MyThemeLine_grid +
      scale_fill_manual(name="Scenario",values=scenariopalette_low2) +
      ggtitle(label=variableplot[i]) +
      guides(color=guide_legend(ncol=1)) +
      ylim(minv-abs(minv)*0.1, maxv+abs(maxv)*0) +
      ylab("") +xlab("") +
      theme(legend.position="none",text=element_text(size=12),  
            axis.text.x=element_text(angle=45, vjust=0.9, hjust=0.5, size = 12)) +
      theme(
        legend.position = "none",
        axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.line = element_line(colour = "white"),
        plot.title = element_text(colour = "white"),
        panel.border = element_blank(),
        axis.ticks=element_blank()
      )
    
    gbox2 <- gbox +facet_wrap(Year~.,ncol=4,strip.position = "bottom") + 
      theme(strip.text.x = element_text(size =10))
    #plot(gbox2)

    outname3 <- paste0("../output/fig/ENGAGE/line_box_FP_3366_cb600/",variableplot[i],".png")  
    multfigure2(gline,gbox2,outname3)
  }
}

# FvsP for 600
for (i in 1:length(variableplot)){
  worlddata3.6 <- filter(worlddata3.5.1,Variable2==variableplot[i], Scenario %in% c("EN_NPi2020_600f","EN_NPi2020_600"))
  worlddata3.6.0 <- filter(worlddata3.6,Year %in% c(2000, 2010, 2020, 2030, 2040, 2050, 2060, 2070, 2080, 2090, 2100))
  worlddata3.6.1 <- ddply(worlddata3.6.0,c("Scenario","Year"),transform,min=quantile(Value,0.33))
  worlddata3.6.2 <- ddply(worlddata3.6.1,c("Scenario","Year"),transform,max=quantile(Value,0.66))
  worlddata3.6.3 <- ddply(worlddata3.6.2,c("Scenario","Year"),transform,v50=quantile(Value,0.50))
  
  
  if(nrow(worlddata3.6)>=2){
    minv <- min(worlddata3.6[7])
    maxv <- max(worlddata3.6[7])
    
    gline <- ggplot(data=worlddata3.6.3) +  
      geom_ribbon(aes(x=Year,ymin=min, ymax=max, fill=Scenario), alpha=0.2) +
      #      geom_line(data=worlddata3.6, aes(x=Year, y = Value, group=interact, colour=interact2), size=0.1,alpha=0.7) + 
      geom_line(aes(x=Year, y = v50, group=interact, colour=Scenario), size=1,alpha=0.7) + 
      #      geom_line(data=worlddata3.6.3, aes(x=Year, y = v50, group=interact, colour=W2100), size=0.1,alpha=0.7) + 
      geom_point(data=worlddata3.6, aes(x=Year, y= Value, colour=Scenario, shape=ModelN),size=2) +
      ylab(Unitlist2$Unit[i]) + xlab("Years") +labs(fill="")+ MyThemeLine_grid2 + 
      ylim(minv-abs(minv)*0.1, maxv+abs(maxv)*0) +
      scale_shape_manual(values=shapepallet) +
      scale_fill_manual(name="Scenario",values=scenariopalette_low3) +
      scale_colour_manual(name="Scenario",values=scenariopalette_low3) + 
      theme(
        legend.position="none", 
        text=element_text(size=14),  
        plot.title = element_text(size=12),
        axis.text.x=element_text(angle=0, vjust=0.9, hjust=0.5, size = 12)) +
      ggtitle(label=variableplot[i]) +
      guides(color=guide_legend(ncol=1))
    plot(gline)
    #    outname <- paste0("../output/fig/ENGAGE/line_box_FP/legend.png")  
    #    ggsave(gline, file=outname, dpi = 600, width=9, height=5,limitsize=FALSE)
    
    worlddata3.7 <- filter(worlddata3.6,Year %in% c(2050,2100))
    gbox <- ggplot() +  
      #      geom_boxplot(data=worlddata3.7,aes(x=W2100, y = Value , fill=W2100), size=0.5,alpha=0.7,range=1) +
      geom_boxplot(data=worlddata3.7,aes(x=Scenario, y = Value , fill=Scenario), size=0.5,alpha=0.7) +
      ylab(Unitlist2$Unit[i]) + xlab("Scenarios") +labs(fill="")+ MyThemeLine_grid +
      scale_fill_manual(name="Scenario",values=scenariopalette_low3) +
      ggtitle(label=variableplot[i]) +
      guides(color=guide_legend(ncol=1)) +
      ylim(minv-abs(minv)*0.1, maxv+abs(maxv)*0) +
      ylab("") +xlab("") +
      theme(legend.position="none",text=element_text(size=12),  
            axis.text.x=element_text(angle=45, vjust=0.9, hjust=0.5, size = 12)) +
      theme(
        legend.position = "none",
        axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        axis.line = element_line(colour = "white"),
        plot.title = element_text(colour = "white"),
        panel.border = element_blank(),
        axis.ticks=element_blank()
      )
    
    gbox2 <- gbox +facet_wrap(Year~.,ncol=4,strip.position = "bottom") + 
      theme(strip.text.x = element_text(size =10))
    plot(gbox2)
    
    outname3 <- paste0("../output/fig/ENGAGE/line_box_FP_3366_cb600_sce/",variableplot[i],".png")  
    multfigure2(gline,gbox2,outname3)
  }
}

#-----box figures

for (i in 1:length(variableplot)){
  worlddata3.6 <- filter(worlddata3.5,Variable2==variableplot[i]) 

  if(nrow(worlddata3.6)>=2){
    worlddata3.7 <- filter(worlddata3.6,Year %in% c(2050,2100),W2100 %in% c('1.5Deg','2Deg'))
    gbox <- ggplot() +  
      geom_boxplot(data=worlddata3.7,aes(x=ModelN, y = Value , fill=ModelN), size=0.5,alpha=0.7,outlier.color = NA) +
      ylab(Unitlist2$Unit[i]) + xlab("Models") +labs(fill="")+ MyThemeLine_grid2 +
      #scale_fill_manual(name="Scenario",values=scenariopalette_low2) +
      theme(legend.position="right", text=element_text(size=12),  
            axis.text.x=element_text(angle=90, vjust=0.9, hjust=1, size = 12)) +
      ggtitle(label=variableplot[i]) +
      guides(color=guide_legend(ncol=1))
    
    gbox2 <- gbox +facet_grid(Year ~ interact2, scales="free_y")
#    plot(gbox2)
    outname2 <- paste0("../output/fig/ENGAGE/box_model/",variableplot[i],".png")
    ggsave(gbox2, file=outname2, dpi = 600, width=9, height=8,limitsize=FALSE)
    
    
  }
}


worlddata3.9 <- filter(worlddata3.0,Variable2 %in% c("Primary_Energy_Biomass_Modern","Land_Cover_Cropland_Energy_Crops")) 
worlddata3.9 <- worlddata3.9[c(-5)]
worlddata3.9 <- worlddata3.9 %>%
  spread(key=Variable2, value=Value, fill=NA)
worlddata3.9 <- na.omit(worlddata3.9)

gschatter2 <- ggplot(worlddata3.9,aes(x=Land_Cover_Cropland_Energy_Crops, y = Primary_Energy_Biomass_Modern, colour=Full_Peak)) + 
  geom_point(aes(shape=ModelN), size=1, alpha=0.7,outlier.color = NA) + 
  xlab("Cropland area for energy crops [million ha]") +
  ylab("Primary energy biomass modern [EJ/year]") + labs(fill="")+ MyThemeLine_grid + 
  scale_colour_manual(name="Budget cap",values=Engagepal) +
  scale_fill_manual(name="Budget cap",values=Engagepal) +
  scale_shape_manual(values=shapepallet) +
  theme(text=element_text(size=12),  
        axis.text.x=element_text(angle=0, vjust=0.9, hjust=1, size = 12)) +
  guides(color=guide_legend(ncol=1))

plot(gschatter2)

outname <- paste0("../output/fig/ENGAGE/schatter/Bio_Land.png")
ggsave(gschatter2, file=outname, dpi = 600, width=9, height=9,limitsize=FALSE)



#--- schatter vs NetZeroYear

for (i in 1:length(variableplot)){
  worlddata3.10 <- filter(Worlddata1.1_CumCO2_NetZero,Variable2==variableplot[i]) 
  varname <- variableplot[i]
  if(nrow(worlddata3.1)>=2){
    worlddata3.11 <- filter(worlddata3.10,Year %in% c(2030,2050,2100))
    
    gschatter3 <- ggplot(worlddata3.11,aes(x=NetZeroYear, y = Value, colour=Year)) +  
      geom_point(aes(shape=ModelN), size=1, alpha=0.7,outlier.color = NA) + 
      ylab(Unitlist2$Unit[i]) + 
      xlab("Year of netzero CO2 emissions [year]") +labs(fill="")+ MyThemeLine_grid + 
#      scale_colour_manual(name="Budget cap",values=Engagepal) +
#      scale_fill_manual(name="Budget cap",values=Engagepal) +
#      scale_shape_manual(values=shapepallet) +
      theme(text=element_text(size=12),  
            axis.text.x=element_text(angle=0, vjust=0.9, hjust=1, size = 12)) +
      ggtitle(label=variableplot[i]) +
      guides(color=guide_legend(ncol=1))
    #gschatter4 <- gschatter3 +facet_wrap(Year~.,ncol=4)
    #plot(gschatter4)
    figurelist[[variableplot[i]]] <- gschatter3 +facet_wrap(Year~.,ncol=4)
    outname <- paste0("../output/fig/ENGAGE/schatter_netzero/",variableplot[i],".png")
    ggsave(figurelist[[variableplot[i]]], file=outname, dpi = 600, width=9, height=4,limitsize=FALSE)

  }
}



#-------AR6 figure
#Cumulative emissions and variables relationship
for (i in 1:length(variableplot)){
  Worlddata1.1_CumCO2.sel <- filter(Worlddata1.1_CumCO2,Year %in% c(2030,2050,2070,2100) & Variable2==variableplot[i]) 
  if(nrow(Worlddata1.1_CumCO2.sel)>=2){
    maxv <- max(Worlddata1.1_CumCO2.sel$Value)
    minv <- min(Worlddata1.1_CumCO2.sel$Value)

    Worlddata1.1_CumCO2.sel <- filter(Worlddata1.1_CumCO2.sel,Year %in% c(2030,2050)) 
    g3 <- ggplot() +
      geom_point(data=Worlddata1.1_CumCO2.sel,aes(x=TCRECumCO2Emis,y=Value,color=Year,fill=Year),shape=21)     + xlim(500,5000)    + 
      ylab(Unitlist2$Unit[i]) + xlab(expression(paste("Cumulative ",CO[2]," emissions until the emissions become zero (Gt) "))) +
      MyThemeLine_grid +labs(size=expression("Carbon Price \n in 2100")) + 
      annotate("rect",xmin=570,xmax=1080,ymin=minv,ymax=maxv,alpha=0.2,fill=spectpal[10]) +
      annotate("rect",xmin=1320,xmax=2270,ymin=minv,ymax=maxv,alpha=0.2,fill=spectpal[8]) +
      ggtitle(label=variableplot[i]) +
      scale_colour_manual(values=pastelpal1)
    outname <- paste0("../output/fig/CumCO2_single/",variableplot[i],".png")
    ggsave(g3, file=outname, dpi = 600, width=5, height=3.5,limitsize=FALSE)
  }
} 

#Periodic figure with AR6 category C1-7
scenariometa <- read.table("../define/meta.txt", header=F, sep="\t")
names(scenariometa) <- c("Scenario","SceMet")
worlddata2.0 <- inner_join(worlddata1.0,scenariometa,by="Scenario")
worlddata2.0 <- worlddata2.0[!grepl("X", worlddata2.0$Year),]
worlddata2.0$Year <- as.numeric(levels(worlddata2.0$Year))[worlddata2.0$Year]
worlddata2.0$SceMet <- factor(worlddata2.0$SceMet, levels = c("no-climate-assessment","C7: above 3.0C","C6: below 3.0C","C5: below 2.5C","C4: higher 2C","C3: lower 2C","C2: 1.5C with high OS","C1: 1.5C with no or low OS"))
worlddata2.0 <- filter(worlddata2.0,SceMet!=c("reference") & SceMet!=c("no-climate-assessment"))

for (i in 1:length(variableplot)){
  worlddata2.1 <- filter(worlddata2.0,Variable2==variableplot[i]) 
  if(nrow(worlddata2.1)>=2){
    worlddata2.2 <- filter(worlddata2.1,Year %in% c(2010,2030,2050,2100))
    gbox <- ggplot() +  
      geom_boxplot(data=worlddata2.2,aes(x=SceMet, y = Value , fill=SceMet), size=0.5,alpha=0.7,outlier.color = NA) + 
      ylab(Unitlist2$Unit[i]) + xlab("Scenarios") +labs(fill="")+ MyThemeLine_grid + 
      scale_fill_manual(values=AR6pal) +
      theme(legend.position="right", text=element_text(size=12),  
            axis.text.x=element_text(angle=45, vjust=0.9, hjust=1, size = 12)) +
      ggtitle(label=variableplot[i]) +
      guides(color=guide_legend(ncol=1))
    gbox2 <- gbox +facet_wrap(Year~.,ncol=4)
    outname <- paste0("../output/fig/box/",variableplot[i],".png")
    ggsave(gbox2, file=outname, dpi = 600, width=9, height=5,limitsize=FALSE)
    
  }
}

#Multiple figures
variableplotM_load <- read.table("../define/variablelistM.txt", header=F, sep=",") 
variableplotM <- as.vector(variableplotM_load$V1)
UnitlistM <- unique(select(worlddata1.0,c("Variable2","Unit")))
UnitlistM2 <- left_join(rename(variableplotM_load,"Variable2"=V1),UnitlistM)

for (i in 1:length(variableplotM)){
  worlddata5.1 <- filter(worlddata2.0,Variable2==variableplotM[i])
  F_name <- paste('F', i, sep = '')
  if(nrow(worlddata5.1)>=2){
    worlddata5.2 <- filter(worlddata5.1,Year %in% c(2010,2030,2050,2100))
    gbox <- ggplot() +  
      geom_boxplot(data=worlddata5.2,aes(x=SceMet, y = Value , fill=SceMet), size=0.5,alpha=0.7,outlier.color = NA) + 
      ylab(UnitlistM2$Unit[i]) + xlab("Scenarios") +labs(fill="")+ MyThemeLine_grid2 + 
      scale_fill_manual(values=AR6pal) +
      theme(legend.position="right", text=element_text(size=9), axis.title.x = element_blank(),   
            element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank()) +
      ggtitle(label=variableplotM[i]) +
      guides(color=guide_legend(ncol=1))
    gbox2 <- gbox + facet_wrap(Year~.,ncol=4)
    YYY <- gtable::gtable_filter(ggplotGrob(gbox2), pattern = "guide-box")
    gbox3 <- gbox2 + theme(legend.position="none")
    assign(F_name, gbox3)
  }
}

XXX <- plot_grid(F1, F2, F3, F4, ncol=2, rel_widths =c(1,1))
MultiFifig <- plot_grid(XXX, YYY, ncol=2, rel_widths=c(5,1))
outname <- paste0("../output/fig/box_multi/MultiFig.png")
ggsave(MultiFifig, file=outname, dpi = 600, width=10, height=5, limitsize=FALSE)


#-----NetZeroYear
NetZeroYear0.1 <- filter(NetZeroYear, Region %in% c("World"))
NetZeroYear1.0 <- ddply(NetZeroYear0.1,c("Scenario"),transform,v50=quantile(NetZeroYear,0.5))
NetZeroYear3.0 <- NetZeroYear1.0 %>% inner_join(scecbfpmap,by="Scenario") %>% inner_join(scenariolist_engage,by=c("Scenario"))

NetZeroYear3.1 <- filter(NetZeroYear3.0,CarbonBudget %in% c(600))
NetZeroYear3.1 <- NetZeroYear3.1[, c(8,10)] 
NetZeroYear3.1$v50 <- as.numeric(NetZeroYear3.1$v50)
NetZeroYear3.2 <- fprename(NetZeroYear3.1)

AFOLUNetZeroYear0.1 <- filter(AFOLUNetZeroYear, Region %in% c("World"))
AFOLUNetZeroYear1.0 <- ddply(AFOLUNetZeroYear0.1,c("Scenario"),transform,v50=quantile(AFOLUNetZeroYear,0.5))
AFOLUNetZeroYear3.0 <- AFOLUNetZeroYear1.0 %>% inner_join(scecbfpmap,by="Scenario") %>% inner_join(scenariolist_engage,by=c("Scenario"))

AFOLUNetZeroYear3.1 <- filter(AFOLUNetZeroYear3.0,CarbonBudget %in% c(600))
AFOLUNetZeroYear3.1 <- AFOLUNetZeroYear3.1[, c(8,10)] 
AFOLUNetZeroYear3.1$v50 <- as.numeric(AFOLUNetZeroYear3.1$v50)
AFOLUNetZeroYear3.2 <- fprename(AFOLUNetZeroYear3.1)

#-----Box for NetZeroYear
AFOLUNetZeroYear0.2 <- AFOLUNetZeroYear
names(AFOLUNetZeroYear0.2) <- c("ModelN","Scenario","Region","Variable3","Unit","NetZeroYear","Value")	
NetZeroYearall<-rbind(NetZeroYear,AFOLUNetZeroYear0.2)

levels(NetZeroYearall$Variable3)[levels(NetZeroYearall$Variable3)=="Emi_CO2"] <- "CO2NetZero"
levels(NetZeroYearall$Variable3)[levels(NetZeroYearall$Variable3)=="Emi_TOT_AFOLU"] <- "AFOLUNetZero"
NetZeroYearall <- NetZeroYearall %>% inner_join(scecbfpmap,by="Scenario") %>% anti_join(scenariomodel_unused,by=c("Scenario","ModelN"))
interact <- interaction(NetZeroYearall$Variable3,NetZeroYearall$Full_Peak, sep=" : ")
NetZeroYearall$Region <- factor(NetZeroYearall$Region, levels=c('REF','MEA','LAM','ASIA','OECD90','World'))
#NetZeroYearall$Region <- factor(NetZeroYearall$Region, levels=c('REF','OAS','SAS','China','MEA','SSA','LAM','POEC','EU28','NAM','World'))
NetZeroYearall <- filter(NetZeroYearall, CarbonBudget %in% c('500','600','700','800','900','1000','1200','1400','1600','1800','2000'))
NetZeroYearall$CarbonBudget <- factor(NetZeroYearall$CarbonBudget, levels=c('300','400','500','600','700','800','900','1000','1200','1400','1600','1800','2000'))


# NetZero year : CO2 vs AFOLU & Full vs peak
gbox <- ggplot() +  
  geom_boxplot(data=NetZeroYearall,aes(x=NetZeroYear, y = Region, fill=interact), size=0.5,alpha=0.7) +
  ylab("") + xlab("Year of Net Zero") +labs(fill="")+ MyThemeLine_grid +
  ggtitle("Regional timing of Net Zero in Full") +
  theme(legend.position="right",text=element_text(size=12),
        legend.text = element_text(size = 10),
        axis.text.x=element_text(angle=0, vjust=0.9, hjust=0.5, size = 13)
        ) 
plot(gbox)
outname <- paste0("../output/fig/ENGAGE/box/AFOLUNetZeroYear_R5.png")  
ggsave(gbox, file=outname, dpi = 600, width=6, height=4, limitsize=FALSE)

# AFOLU NetZero year : Region, Full vs peak

AFOLUNetZeroYear4.0 <- filter(NetZeroYearall, Variable3 %in% c("AFOLUNetZero"))
#AFOLUNetZeroYear4.1 <- filter(AFOLUNetZeroYear4.0, CarbonBudget %in% c("600"))

gboxafolur <- ggplot() +  
  geom_boxplot(data=AFOLUNetZeroYear4.0,aes(x=NetZeroYear, y = Region, fill=Full_Peak), size=0.5,alpha=0.7) +
  ylab("") + xlab("Year of Net Zero") +labs(fill="") + MyThemeLine_grid2h +
  ggtitle("Net Zero year of AFOLU GHG emissions") +
  theme(legend.position="top",text=element_text(size=12),
        legend.text = element_text(size = 10),
        axis.text.x=element_text(angle=0, vjust=0.9, hjust=0.5, size = 13)
  )  +
  scale_fill_manual(values=pastelpal1)

plot(gboxafolur)
outname <- paste0("../output/fig/ENGAGE/box/AFOLUNetZeroYear_FullvsPeak_region.pdf")  
ggsave(gboxafolur, file=outname, dpi = 600, width=4.5, height=4, limitsize=FALSE)
outname <- paste0("../output/fig/ENGAGE/pdf/Fig2c.pdf")
ggsave(gboxafolur, file=outname, dpi = 600, width=4.5, height=4, limitsize=FALSE)

# AFOLU NetZero year : by CB, Full vs peak

AFOLUNetZeroYear4.2 <- filter(AFOLUNetZeroYear4.0, Region %in% c("World"))

gboxafolu <- ggplot() +  
  geom_boxplot(data=AFOLUNetZeroYear4.2,aes(x=NetZeroYear, y = CarbonBudget, fill=Full_Peak), size=0.5,alpha=0.7) +
  ylab("") + xlab("Year of Net Zero") +labs(fill="") + MyThemeLine_grid2h +
  ggtitle("Net Zero year of AFOLU GHG emissions") +
  theme(legend.position="top",text=element_text(size=12),
        legend.text = element_text(size = 10),
        axis.text.x=element_text(angle=0, vjust=0.9, hjust=0.5, size = 13)
  ) +
  scale_fill_manual(values=pastelpal1)
plot(gboxafolu)
outname <- paste0("../output/fig/ENGAGE/box/AFOLUNetZeroYear_FullvsPeak_cb.png")  
ggsave(gboxafolu, file=outname, dpi = 600, width=4.5, height=4, limitsize=FALSE)


# by Region, Peak only
#NetZeroYearallpeak <- filter(NetZeroYearall, Full_Peak %in% c('Peak') & CarbonBudget %in% c('600'))
NetZeroYearallpeak <- filter(NetZeroYearall, Full_Peak %in% c('Peak'))

gboxpeak <- ggplot() +  
  geom_boxplot(data=NetZeroYearallpeak,aes(x=NetZeroYear, y = Region, fill=Variable3), size=0.5,alpha=0.7) +
  ylab("") + xlab("Year of Net Zero") +labs(fill="")+ MyThemeLine_grid +
  ggtitle("Regional timing of Net Zero in Peak") +
  theme(legend.position="top",text=element_text(size=12),
        legend.text = element_text(size = 10),
        axis.text.x=element_text(angle=0, vjust=0.9, hjust=0.5, size = 13)
  ) 

plot(gboxpeak)

outname <- paste0("../output/fig/ENGAGE/box/AFOLUNetZeroYear_byRegion_Peakonly_R5rev.pdf")  
ggsave(gboxpeak, file=outname, dpi = 600, width=4.5, height=4, limitsize=FALSE)
outname <- paste0("../output/fig/ENGAGE/pdf/Fig2b.pdf")  
ggsave(gboxpeak, file=outname, dpi = 600, width=4.5, height=4, limitsize=FALSE)

# NetZeroYear vs CarbonBudget
NetZeroYearallpeakworld <- filter(NetZeroYearallpeak, Region %in% c('World') & CarbonBudget %in% c('500','600','700','800','900','1000','1200','1400','1600','1800','2000'))
NetZeroYearallpeakworld$CarbonBudget <- factor(NetZeroYearallpeakworld$CarbonBudget, levels=c('300','400','500','600','700','800','900','1000','1200','1400','1600','1800','2000'))

gboxfullcb <- ggplot() +  
  geom_boxplot(data=NetZeroYearallpeakworld,aes(x=NetZeroYear, y = CarbonBudget, fill=Variable3), size=0.5,alpha=0.7) +
  ylab("Carbon budget [GtCO2]") + xlab("Year of Net Zero") +labs(fill="")+ MyThemeLine_grid +
  ggtitle("Global timing of Net Zero in Peak") +
  theme(legend.position="top",text=element_text(size=12),
        legend.text = element_text(size = 10),
        axis.text.x=element_text(angle=0, vjust=0.9, hjust=0.5, size = 13)
  ) 

plot(gboxfullcb)

outname <- paste0("../output/fig/ENGAGE/box/AFOLUNetZeroYear_Peakonly_R5.pdf")  
ggsave(gboxfullcb, file=outname, dpi = 600, width=4.5, height=4, limitsize=FALSE)
outname <- paste0("../output/fig/ENGAGE/pdf/Fig2a.pdf")  
ggsave(gboxfullcb, file=outname, dpi = 600, width=4.5, height=4, limitsize=FALSE)


#-----

fvarname1 <- function(y){
  levels(y$Variable2)[levels(y$Variable2)=="Land_Cover_Pasture_Change_2010"] <- "Pasture"
  levels(y$Variable2)[levels(y$Variable2)=="Land_Cover_Forest_Change_2010"] <- "Forest"
  levels(y$Variable2)[levels(y$Variable2)=="Land_Cover_Cropland_NonEnergy_Crops_Change_2010"] <- "Cropland for NonEnergy"
  levels(y$Variable2)[levels(y$Variable2)=="Land_Cover_Cropland_Energy_Crops_Change_2010"] <- "Cropland for Energy"
  levels(y$Variable2)[levels(y$Variable2)=="Land_Cover_Other_Natural_Land_Change_2010"] <- "Other_Natural_Land"
  levels(y$Variable2)[levels(y$Variable2)=="Land_Pos_Change"] <- "Land_Pos_Change"
  levels(y$Variable2)[levels(y$Variable2)=="Land_Neg_Change"] <- "Land_Neg_Change"
  
  levels(y$Variable2)[levels(y$Variable2)=="Carbon_Sequestration_CCS_Biomass"] <- "Carbon_Sequestration_BECCS"
  levels(y$Variable2)[levels(y$Variable2)=="Emi_CO2_AFOLU"] <- "CO2_AFOLU"
  levels(y$Variable2)[levels(y$Variable2)=="Emi_CH4_AFOLU"] <- "CH4_AFOLU"
  levels(y$Variable2)[levels(y$Variable2)=="Emi_N2O_AFOLU"] <- "N2O_AFOLU"
  levels(y$Variable2)[levels(y$Variable2)=="Emi_TOT_AFOLU"] <- "NET_AFOLU"
  levels(y$Variable2)[levels(y$Variable2)=="Emi_TOT_AFOLU_NoBECCS"] <- "NET_AFOLU_NoBECCS"
  levels(y$Variable2)[levels(y$Variable2)=="Emi_Pos_AFOLU"] <- "Pos_AFOLU"
  levels(y$Variable2)[levels(y$Variable2)=="Emi_Neg_AFOLU"] <- "Neg_AFOLU"
  levels(y$Variable2)[levels(y$Variable2)=="Emissions_CH4_AFOLU_Agriculture_Livestock_Manure_Management"] <- "CH4_Manure_Management"
  levels(y$Variable2)[levels(y$Variable2)=="Emissions_CH4_AFOLU_Agriculture_Livestock_Enteric_Fermentation"] <- "CH4_Enteric_Fermentation"
  levels(y$Variable2)[levels(y$Variable2)=="Emissions_CH4_AFOLU_Agriculture_Rice"] <- "CH4_Rice"
  levels(y$Variable2)[levels(y$Variable2)=="Emissions_N2O_AFOLU_Agriculture_Livestock_Manure_Management"] <- "N2O_Manure_Management"
  levels(y$Variable2)[levels(y$Variable2)=="Emissions_N2O_AFOLU_Agriculture_Managed_Soils"] <- "N2O_Managed_Soils"

  levels(y$Variable2)[levels(y$Variable2)=="Price_Agriculture_Non-Energy_Cops_and_Livestock_Index"] <- "Agricultural\nprice\n[2005=1]"
  levels(y$Variable2)[levels(y$Variable2)=="Food_Demand"] <- "Food\nconsumption\n[kcal/person/day]"
  levels(y$Variable2)[levels(y$Variable2)=="Food_Demand_Livestock"] <- "Livestock\nproducts\nconsumption\n[kcal/person/day]"
  levels(y$Variable2)[levels(y$Variable2)=="Population_Risk_of_Hunger"] <- "Popuation\nat risk\nof hunger\n[million\npeople]"
  levels(y$Variable2)[levels(y$Variable2)=="Yield_Cereal"] <- "Yield\ncereal\n[t DM/ha/yr]"
  levels(y$Variable2)[levels(y$Variable2)=="Primary_Energy_Biomass_Modern"] <- "Primary\nbioenergy\n[EJ/yr]"
  levels(y$Variable2)[levels(y$Variable2)=="Price_Carbon"] <- "Carbon\nprice\n[US$2010/t CO2]"
  levels(y$Variable2)[levels(y$Variable2)=="Fertilizer_Use_Nitrogen"] <- "Nitrogen\nfertilizer\nuse\n[Tg N/yr]"
  levels(y$Variable2)[levels(y$Variable2)=="Water_Withdrawal_Irrigatio"] <- "Water withdrawal\nfor irrigation\n[km3/year]"
                                            
  z <- y
  return(z)
}

fprename <- function(y){
  levels(y$Full_Peak)[levels(y$Full_Peak)=="Full"] <- "EOC"
  levels(y$Full_Peak)[levels(y$Full_Peak)=="Peak"] <- "NZ"

  z <- y
  return(z)
}

modelrename <- function(y){
  levels(y$ModelN)[levels(y$ModelN)=="AIM_CGE 2_2"] <- "AIM/CGE"
  levels(y$ModelN)[levels(y$ModelN)=="COFFEE 1_1"] <- "COFFEE"
  levels(y$ModelN)[levels(y$ModelN)=="IMAGE 3_0"] <- "IMAGE"
  levels(y$ModelN)[levels(y$ModelN)=="MESSAGEix-GLOBIOM_1.0"] <- "MESSAGEix-GLOBIOM"
  levels(y$ModelN)[levels(y$ModelN)=="POLES_ENGAGE"] <- "POLES"
  levels(y$ModelN)[levels(y$ModelN)=="REMIND-MAgPIE 2_0-4_1"] <- "REMIND-MAgPIE"
  levels(y$ModelN)[levels(y$ModelN)=="WITCH 5_0"] <- "WITCH"
  
  
  z <- y
  return(z)
}

#-----area by Emissions sources

AFOLUlist_load <- read.table("../define/AFOLUlist.txt", header=F, sep=",") 
AFOLUlist <- as.vector(AFOLUlist_load$V1) 
AFOLUlist2 <- read.table("../define/AFOLUlist2.txt", header=T, sep=",") 
names(AFOLUlist2) <- c("Variable2")
AFOLUlist3 <- read.table("../define/AFOLUlist3.txt", header=T, sep=",") 
names(AFOLUlist3) <- c("Variable2")
AFOLUlist4 <- read.table("../define/AFOLUlist4.txt", header=T, sep=",") 
names(AFOLUlist4) <- c("Variable2")


#worlddata_AFOLU1.3$Year <- as.numeric(worlddata_AFOLU1.3$Year)


worlddata_AFOLU <- filter(worlddata3.0,Region %in% c("World"),Variable2 %in% AFOLUlist,Year %in% c(2000, 2010, 2020, 2030, 2040, 2050, 2060, 2070, 2080, 2090, 2100))
worlddata_AFOLU$Value <- ifelse(worlddata_AFOLU$Unit=="Mt CH4/yr", worlddata_AFOLU$Value*25, worlddata_AFOLU$Value) 
worlddata_AFOLU$Value <- ifelse(worlddata_AFOLU$Unit=="kt N2O/yr", worlddata_AFOLU$Value*298/1000, worlddata_AFOLU$Value) 
worlddata_AFOLU$Unit <- "Mt CO2eq/yr"
worlddata_AFOLU$Full_Peak <- factor(worlddata_AFOLU$Full_Peak, levels=c('Peak','Full'))

worlddata_AFOLU0.1 <- filter(worlddata_AFOLU,Year %in% c(2050)) %>% anti_join(AFOLUlist4,by=c("Variable2"))
worlddata_AFOLU0.2 <- filter(worlddata_AFOLU,CarbonBudget %in% c(600)) %>% anti_join(AFOLUlist4,by=c("Variable2"))
worlddata_AFOLU0.3 <- filter(worlddata_AFOLU0.1,CarbonBudget %in% c(600),Full_Peak %in% c('Peak'))
worlddata_AFOLU0.4 <- filter(worlddata_AFOLU,CarbonBudget %in% c(600)) %>% anti_join(AFOLUlist4,by=c("Variable2"))
#worlddata_AFOLU0.5 <- filter(worlddata_AFOLU0.2,ModelN %in% c('AIM_CGE 2_2'))

worlddata_AFOLU1.0 <- ddply(worlddata_AFOLU,c("Variable2","Full_Peak","CarbonBudget","Year"),transform,v50=quantile(Value,0.5))
worlddata_AFOLU1.0a <- ddply(worlddata_AFOLU1.0,c("Variable2","Full_Peak","CarbonBudget","Year"),transform,min=quantile(Value,0))
worlddata_AFOLU1.0a <- ddply(worlddata_AFOLU1.0a,c("Variable2","Full_Peak","CarbonBudget","Year"),transform,max=quantile(Value,1))

worlddata_AFOLU1.0.1 <- worlddata_AFOLU1.0a[worlddata_AFOLU1.0a$ModelN=="MESSAGEix-GLOBIOM_1.0", c(1,2,3,4,5,6,7,8,9,10,11,12)] 
worlddata_AFOLU1.0.1 <- fvarname1(worlddata_AFOLU1.0.1)
worlddata_AFOLU1.0.1 <- fprename(worlddata_AFOLU1.0.1)

worlddata_AFOLU1.1a <- worlddata_AFOLU1.0.1 %>% filter(Year %in% c(2050)) %>% anti_join(AFOLUlist2,by=c("Variable2"))
worlddata_AFOLU1.1b <- worlddata_AFOLU1.0.1 %>% filter(Year %in% c(2050)) %>% inner_join(AFOLUlist3,by=c("Variable2"))
worlddata_AFOLU1.2a <- worlddata_AFOLU1.0.1 %>% filter(CarbonBudget %in% c(600)) %>% anti_join(AFOLUlist2,by=c("Variable2"))
worlddata_AFOLU1.2b <- worlddata_AFOLU1.0.1 %>% filter(CarbonBudget %in% c(600)) %>% inner_join(AFOLUlist3,by=c("Variable2"))

worlddata_AFOLU1.2a2<-worlddata_AFOLU1.2a
worlddata_AFOLU1.2a$Variable2 <- factor(worlddata_AFOLU1.2a$Variable2, levels=c('N2O_AFOLU','CH4_AFOLU','CO2_AFOLU','Carbon_Sequestration_BECCS'))
worlddata_AFOLU1.2a2$Variable2 <- factor(worlddata_AFOLU1.2a2$Variable2, levels=c('N2O_AFOLU','CH4_AFOLU','Carbon_Sequestration_BECCS','CO2_AFOLU'))

worlddata_AFOLU1.3 <- worlddata_AFOLU1.2a[c(-7,-11,-12)]
worlddata_AFOLU1.3.1 <- worlddata_AFOLU1.3 %>% spread(key=Variable2, value=v50, fill=0)

worlddata_AFOLU1.3.1[12] <- worlddata_AFOLU1.3.1[8]+worlddata_AFOLU1.3.1[9]+worlddata_AFOLU1.3.1[10]+worlddata_AFOLU1.3.1[11]
worlddata_AFOLU1.3.1[13] <- worlddata_AFOLU1.3.1[8]+worlddata_AFOLU1.3.1[9]+worlddata_AFOLU1.3.1[10]
names(worlddata_AFOLU1.3.1) <- c("ModelN","Scenario","Region","Unit","Year","CarbonBudget","Full_Peak",'N2O_AFOLU','CH4_AFOLU','CO2_AFOLU','Carbon_Sequestration_BECCS','NET_AFOLU','NET_AFOLU_NoBECCS')	#Rename the column
worlddata_AFOLU1.3.1 <- worlddata_AFOLU1.3.1 %>% gather(key=Variable2,value=v50,-ModelN,-Scenario,-Region,-Unit,-Year,-CarbonBudget,-Full_Peak)
worlddata_AFOLU1.4 <- filter(worlddata_AFOLU1.3.1,Variable2==c('NET_AFOLU'))
worlddata_AFOLU1.5 <- filter(worlddata_AFOLU1.3.1,Variable2==c('NET_AFOLU_NoBECCS'))

g10 <- ggplot() +
  geom_bar(data=subset(worlddata_AFOLU1.1a,v50>0), aes(x=CarbonBudget,y=v50/1000,fill=Variable2), position="stack", stat="identity") +
  geom_bar(data=subset(worlddata_AFOLU1.1a,v50<0), aes(x=CarbonBudget,y=v50/1000,fill=Variable2), position="stack", stat="identity") +
  geom_errorbar(data=subset(worlddata_AFOLU1.1b,v50>0), aes(x=CarbonBudget,ymin=min/1000, ymax=max/1000), colour='gray80', width=.2,position=position_dodge(.9)) + 
  geom_errorbar(data=subset(worlddata_AFOLU1.1b,v50<0), aes(x=CarbonBudget,ymin=min/1000, ymax=max/1000), colour='gray80', width=.2,position=position_dodge(.9)) + 
  geom_hline(yintercept=0,linetype="dashed") +
  ylab("Gt CO2eq/yr") + xlab("Carbon budget [GtCO2]") + MyThemeLine_grid2h + 
  ggtitle(label="GHG Emission and Sequestration in AFOLU in 2050") +
  theme(legend.position="right", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10)) +
  scale_fill_manual(values=pastelpal1)
g10.1 <- g10+facet_wrap(Full_Peak~.)
plot(g10.1)
outname <- paste0("../output/fig/ENGAGE/area/Emission_AFOLU_bysource_cb.pdf")
ggsave(g10.1, file=outname, dpi = 600, width=8, height=4,limitsize=FALSE)
outname <- paste0("../output/fig/ENGAGE/pdf/Fig1e.pdf")
ggsave(g10.1, file=outname, dpi = 600, width=8, height=4,limitsize=FALSE)

# by models
g11 <- ggplot() +
  geom_bar(data=subset(worlddata_AFOLU0.1,Value>0), aes(x=CarbonBudget,y=Value/1000,fill=Variable2), position="stack", stat="identity") +
  geom_bar(data=subset(worlddata_AFOLU0.1,Value<0), aes(x=CarbonBudget,y=Value/1000,fill=Variable2), position="stack", stat="identity") +
  geom_hline(yintercept=0,linetype="dashed") +
  ylab("Gt CO2eq/yr") + xlab("Carbon budget [GtCO2]") + MyThemeLine_grid2 +
  ggtitle(label="GHG Emission and Sequestration in AFOLU in 2050") +
  theme(legend.position="right", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10)) +
  scale_fill_manual(values=pastelpal1)
g11.1 <- g11+facet_grid(ModelN~Full_Peak)
plot(g11.1)
outname <- paste0("../output/fig/ENGAGE/area/Emission_AFOLU_bysource_model.png")
ggsave(g11.1, file=outname, dpi = 600, width=8, height=8,limitsize=FALSE)

# by Time series
g12 <- ggplot() +
#  geom_bar(data=subset(worlddata_AFOLU1.2a,v50>0), aes(x=Year,y=v50/1000,fill=Variable2), position="stack", stat="identity") +
#  geom_bar(data=subset(worlddata_AFOLU1.2a2,v50<0), aes(x=Year,y=v50/1000,fill=Variable2), position="stack", stat="identity") +
  geom_area(data=subset(worlddata_AFOLU1.2a,v50>0), aes(x=Year,y=v50/1000,fill=Variable2)) +
  geom_area(data=subset(worlddata_AFOLU1.2a2,v50<0), aes(x=Year,y=v50/1000,fill=Variable2)) +
  geom_errorbar(data=subset(worlddata_AFOLU1.2b,v50>0), aes(x=Year,ymin=min/1000, ymax=max/1000), colour='gray80', width=.2,position=position_dodge(.9)) + 
  geom_errorbar(data=subset(worlddata_AFOLU1.2b,v50<0), aes(x=Year,ymin=min/1000, ymax=max/1000), colour='gray80', width=.2,position=position_dodge(.9)) + 
  geom_line(data=worlddata_AFOLU1.4, aes(x=Year,y=v50/1000), colour="black") +
  geom_line(data=worlddata_AFOLU1.5, aes(x=Year,y=v50/1000), colour="black", linetype="dotted") +
  geom_hline(yintercept=0,linetype="solid") +
  geom_vline(data=NetZeroYear3.2,aes(xintercept=v50),colour="blue") +
  geom_vline(data=AFOLUNetZeroYear3.2,aes(xintercept=v50),colour="red") +
  annotate("text", x=2062, y=5, label= "Year of CO2 net zero", color="blue",size=3,angle=90) +
  annotate("text", x=2050, y=5, label= "Year of AFOLU net zero", color="red",size=3,angle=90) +
  ylab("Gt CO2eq/yr") + xlab("Year") + MyThemeLine_grid2h +
  ggtitle(label="GHG Emissions and Sequestration in AFOLU with carbon budget of 600GtCO2") +
  theme(legend.position="right", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10)) +
  scale_fill_manual(values=pastelpal1)
g12.1 <- g12+facet_wrap(Full_Peak~.)
plot(g12.1)
outname <- paste0("../output/fig/ENGAGE/area/Emission_AFOLU_timeseriesv2.pdf")
ggsave(g12.1, file=outname, dpi = 600, width=8, height=4,limitsize=FALSE)
outname <- paste0("../output/fig/ENGAGE/pdf/Fig1a.pdf")
ggsave(g12.1, file=outname, dpi = 600, width=8, height=4,limitsize=FALSE)

# by models
g13 <- ggplot() +
  geom_area(data=subset(worlddata_AFOLU0.2,Value>0), aes(x=Year,y=Value/1000,fill=Variable2)) +
  geom_area(data=subset(worlddata_AFOLU0.2,Value<0), aes(x=Year,y=Value/1000,fill=Variable2)) +
  geom_hline(yintercept=0,linetype="dashed") +
  ylab("Gt CO2eq/yr") + xlab("Year") + MyThemeLine_grid2 +
  ggtitle(label="GHG Emission and Sequestration in AFOLU with carbon budget of 600GtCO2") +
  theme(legend.position="right", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10)) +
  scale_fill_manual(values=pastelpal1)
g13.1 <- g13+facet_grid(ModelN~Full_Peak)
plot(g13.1)
outname <- paste0("../output/fig/ENGAGE/area/Emission_AFOLU_timeseries_model.png")
ggsave(g13.1, file=outname, dpi = 600, width=8, height=8,limitsize=FALSE)

worlddata_AFOLU0.3 <- filter(worlddata_AFOLU0.3, ModelN!=c('GEM-E3_092019'), Variable2!=c('Emi_TOT_AFOLU_NoBECCS'))
worlddata_AFOLU0.3 <- modelrename(worlddata_AFOLU0.3)

# 2050, 600 by models
g14 <- ggplot() +
  geom_bar(data=subset(worlddata_AFOLU0.3,Value>0), aes(x=ModelN,y=Value/1000,fill=Variable2), position="stack", stat="identity") +
  geom_bar(data=subset(worlddata_AFOLU0.3,Value<0), aes(x=ModelN,y=Value/1000,fill=Variable2), position="stack", stat="identity") +
  geom_hline(yintercept=0,linetype="dashed") +
  ylab("Gt CO2eq/yr") + xlab("") + MyThemeLine_grid2 +
  ggtitle(label="AFOLU GHG Emissions for 2050 and 600GtCO2") +
  theme(legend.position="right", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10)) +
  scale_fill_manual(values=pastelpal1)
#g14.1 <- g14+facet_wrap(Full_Peak~.)
plot(g14)
outname <- paste0("../output/fig/ENGAGE/area/Emission_AFOLU_bysource_model_2050_Peakonly.pdf")
ggsave(g14, file=outname, dpi = 600, width=6, height=5,limitsize=FALSE)
outname <- paste0("../output/fig/ENGAGE/pdf/Fig1d.pdf")
ggsave(g14, file=outname, dpi = 600, width=6, height=5,limitsize=FALSE)


# NetEmissions by Time series
g15 <- ggplot() +
  geom_line(data=worlddata_AFOLU0.4, aes(x=Year,y=Value/1000, group=ModelN,colour=ModelN)) +
  geom_hline(yintercept=0,linetype="dashed") +
  geom_vline(data=NetZeroYear3.1,aes(xintercept=v50),colour="blue") +
  geom_vline(data=AFOLUNetZeroYear3.1,aes(xintercept=v50),colour="red") +
  annotate("text", x=2062, y=5, label= "Year of CO2 net zero", color="blue",size=3,angle=90) +
  annotate("text", x=2050, y=5, label= "Year of AFOLU net zero", color="red",size=3,angle=90) +
  ylab("Gt CO2eq/yr") + xlab("Year") + MyThemeLine_grid2 +
  ggtitle(label="GHG Emissions and Sequestration in AFOLU with carbon budget of 600GtCO2") +
  theme(legend.position="right", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10)) +
  scale_fill_manual(values=pastelpal1)
g15.1 <- g15+facet_wrap(Full_Peak~.)
plot(g15.1)
outname <- paste0("../output/fig/ENGAGE/area/AFOLUNetEmissions_timeseries_model.png")
ggsave(g15.1, file=outname, dpi = 600, width=8, height=4,limitsize=FALSE)


#--Area by Land

Landcoverlist_load <- read.table("../define/Landcoverlist_change.txt", header=F, sep=",") 
Landcoverlist <- as.vector(Landcoverlist_load$V1) 
Landcoverlist2 <- read.table("../define/Landcoverlist2.txt", header=T, sep=",") 
names(Landcoverlist2) <- c("Variable2")

worlddata_Land <- filter(worlddata3.0,Region %in% c('World'),Variable2 %in% Landcoverlist,Year %in% c(2000, 2010, 2020, 2030, 2040, 2050, 2060, 2070, 2080, 2090, 2100))
worlddata_Land <- fvarname1(worlddata_Land)
worlddata_Land$Full_Peak <- factor(worlddata_Land$Full_Peak, levels=c('Peak','Full'))
worlddata_Land <- fprename(worlddata_Land)

worlddata_Land0.1 <- filter(worlddata_Land,Year %in% c(2050)) %>% anti_join(Landcoverlist2,by=c("Variable2"))
worlddata_Land0.2 <- filter(worlddata_Land,CarbonBudget %in% c(600)) %>% anti_join(Landcoverlist2,by=c("Variable2"))
worlddata_Land0.3 <- filter(worlddata_Land0.1,CarbonBudget %in% c(600),Full_Peak %in% c('NZ'))

worlddata_Land1.0 <- ddply(worlddata_Land,c("Variable2","Full_Peak","CarbonBudget","Year"),transform,v50=quantile(Value,0.5))
worlddata_Land1.0a <- ddply(worlddata_Land1.0,c("Variable2","Full_Peak","CarbonBudget","Year"),transform,min=quantile(Value,0))
worlddata_Land1.0a <- ddply(worlddata_Land1.0a,c("Variable2","Full_Peak","CarbonBudget","Year"),transform,max=quantile(Value,1))

worlddata_Land1.0.1 <- worlddata_Land1.0a[worlddata_Land1.0a$ModelN=="MESSAGEix-GLOBIOM_1.0", c(-7)]
#worlddata_Land1.0.1 <- fprename(worlddata_Land1.0.1)
worlddata_Land1.1 <- filter(worlddata_Land1.0.1,Year %in% c(2050))
worlddata_Land1.2 <- filter(worlddata_Land1.0.1,CarbonBudget %in% c(600))

worlddata_Land1.1a <- worlddata_Land1.0.1 %>% filter(Year %in% c(2050)) %>% anti_join(Landcoverlist2,by=c("Variable2"))
worlddata_Land1.1b <- worlddata_Land1.0.1 %>% filter(Year %in% c(2050)) %>% inner_join(Landcoverlist2,by=c("Variable2"))
worlddata_Land1.2a <- worlddata_Land1.0.1 %>% filter(CarbonBudget %in% c(600)) %>% anti_join(Landcoverlist2,by=c("Variable2"))
worlddata_Land1.2b <- worlddata_Land1.0.1 %>% filter(CarbonBudget %in% c(600)) %>% inner_join(Landcoverlist2,by=c("Variable2"))


g5 <- ggplot() +
  geom_bar(data=subset(worlddata_Land1.1a,v50>0), aes(x=CarbonBudget,y=v50,fill=Variable2), position="stack", stat="identity") +
  geom_bar(data=subset(worlddata_Land1.1a,v50<0), aes(x=CarbonBudget,y=v50,fill=Variable2), position="stack", stat="identity") +
  geom_errorbar(data=subset(worlddata_Land1.1b,v50>0), aes(x=CarbonBudget,ymin=min, ymax=max), colour='gray80', width=.2,position=position_dodge(.9)) + 
  geom_errorbar(data=subset(worlddata_Land1.1b,v50<0), aes(x=CarbonBudget,ymin=min, ymax=max), colour='gray80', width=.2,position=position_dodge(.9)) + 
  geom_hline(yintercept=0,linetype="dashed") +
  ylab(expression(paste("million ha"))) + xlab("Carbon budget [GtCO2]") + MyThemeLine_grid2h +  
  ggtitle(label="Land Use Change in 2050 with respect to 2010") +
  theme(legend.position="right", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10)) +
  scale_fill_manual(values=fillpalland)
g5.1 <- g5+facet_wrap(Full_Peak~.)
plot(g5.1)
outname <- paste0("../output/fig/ENGAGE/area/LUCwrt2010_2050_cb.pdf")
ggsave(g5.1, file=outname, dpi = 600, width=8, height=4,limitsize=FALSE)
outname <- paste0("../output/fig/ENGAGE/pdf/Fig3e.pdf")
ggsave(g5.1, file=outname, dpi = 600, width=8, height=4,limitsize=FALSE)

g6 <- ggplot() +
  geom_bar(data=subset(worlddata_Land0.1,Value>0), aes(x=CarbonBudget,y=Value,fill=Variable2), position="stack", stat="identity") +
  geom_bar(data=subset(worlddata_Land0.1,Value<0), aes(x=CarbonBudget,y=Value,fill=Variable2), position="stack", stat="identity") +
  geom_hline(yintercept=0,linetype="dashed") +
  ylab(expression(paste("million ha"))) + xlab("Carbon budget [GtCO2]") + MyThemeLine_grid2h +
  ggtitle(label="Land Use Change in 2050 with respect to 2010") +
  theme(legend.position="right", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10)) +
  scale_fill_manual(values=fillpalland)
g6.1 <- g6+facet_grid(ModelN~Full_Peak)
plot(g6.1)
outname <- paste0("../output/fig/ENGAGE/area/LUCwrt2010_2050_cb_model.png")
ggsave(g6.1, file=outname, dpi = 600, width=8, height=8,limitsize=FALSE)

# by Time series
g7 <- ggplot() +
  geom_area(data=subset(worlddata_Land1.2a,v50>0), aes(x=Year,y=v50,fill=Variable2)) +
  geom_area(data=subset(worlddata_Land1.2a,v50<0), aes(x=Year,y=v50,fill=Variable2)) +
  geom_errorbar(data=subset(worlddata_Land1.2b,v50>0), aes(x=Year,ymin=min, ymax=max), colour='gray80', width=.2,position=position_dodge(.9)) + 
  geom_errorbar(data=subset(worlddata_Land1.2b,v50<0), aes(x=Year,ymin=min, ymax=max), colour='gray80', width=.2,position=position_dodge(.9)) + 
  geom_hline(yintercept=0,linetype="dashed") +
  geom_vline(data=NetZeroYear3.2,aes(xintercept=v50),colour="blue") +
  geom_vline(data=AFOLUNetZeroYear3.2,aes(xintercept=v50),colour="red") +
  annotate("text", x=2062, y=5, label= "Year of CO2 net zero", color="blue",size=3,angle=90) +
  annotate("text", x=2050, y=5, label= "Year of AFOLU net zero", color="red",size=3,angle=90) +
  ylab(expression(paste("million ha"))) + xlab("Year") + MyThemeLine_grid2h +
  ggtitle(label="Land Use Change with respect to 2010 with carbon budget of 600GtCO2") +
  theme(legend.position="right", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10)) +
  scale_fill_manual(values=fillpalland)
g7.1 <- g7+facet_wrap(Full_Peak~.)
plot(g7.1)
outname <- paste0("../output/fig/ENGAGE/area/LUC_timeseries.pdf")
ggsave(g7.1, file=outname, dpi = 600, width=8, height=4,limitsize=FALSE)
outname <- paste0("../output/fig/ENGAGE/pdf/Fig3a.pdf")
ggsave(g7.1, file=outname, dpi = 600, width=8, height=4,limitsize=FALSE)

# by models
g8 <- ggplot() +
  geom_area(data=subset(worlddata_Land0.2,Value>0), aes(x=Year,y=Value,fill=Variable2)) +
  geom_area(data=subset(worlddata_Land0.2,Value<0), aes(x=Year,y=Value,fill=Variable2)) +
  geom_hline(yintercept=0,linetype="dashed") +
  ylab(expression(paste("million ha"))) + xlab("Year") + MyThemeLine_grid2h +
  ggtitle(label="Land Use Change with respect to 2010 with carbon budget of 600GtCO2") +
  theme(legend.position="right", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10)) +
  scale_fill_manual(values=fillpalland)
g8.1 <- g8+facet_grid(ModelN~Full_Peak)
plot(g8.1)
outname <- paste0("../output/fig/ENGAGE/area/LUC_timeseries_model.png")
ggsave(g8.1, file=outname, dpi = 600, width=8, height=8,limitsize=FALSE)

worlddata_Land0.3 <- modelrename(worlddata_Land0.3)
# Land for 2050, 600 by models
g20 <- ggplot() +
  geom_bar(data=subset(worlddata_Land0.3,Value>0), aes(x=ModelN,y=Value,fill=Variable2), position="stack", stat="identity") +
  geom_bar(data=subset(worlddata_Land0.3,Value<0), aes(x=ModelN,y=Value,fill=Variable2), position="stack", stat="identity") +
  geom_hline(yintercept=0,linetype="dashed") +
  ylab(expression(paste("million ha"))) + xlab("") + MyThemeLine_grid2 +
  ggtitle(label="Land Use Change in 2050 with respect to 2010 for 600GtCO2") +
  theme(legend.position="right", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10)) +
  scale_fill_manual(values=fillpalland)
plot(g20)
outname <- paste0("../output/fig/ENGAGE/area/LUCwrt2010_2050_model_Peakonly.pdf")
ggsave(g20, file=outname, dpi = 600, width=5, height=5,limitsize=FALSE)
outname <- paste0("../output/fig/ENGAGE/pdf/Fig3d.pdf")
ggsave(g20, file=outname, dpi = 600, width=5, height=5,limitsize=FALSE)



#--Land by Region

regionaldata_Land <- filter(worlddata3.0,Variable2 %in% Landcoverlist)
regionaldata_Land <- filter(regionaldata_Land, Year %in% c(2050), CarbonBudget %in% c(600), Region!=c("World"))
regionaldata_Land <- fvarname1(regionaldata_Land)
regionaldata_Land$Full_Peak <- factor(regionaldata_Land$Full_Peak, levels=c('Peak','Full'))

regionaldata_Land1.0 <- ddply(regionaldata_Land,c("Variable2","Full_Peak","Region"),transform,v50=quantile(Value,0.5))
regionaldata_Land1.0a <- ddply(regionaldata_Land1.0,c("Variable2","Full_Peak","Region"),transform,min=quantile(Value,0))
regionaldata_Land1.0a <- ddply(regionaldata_Land1.0a,c("Variable2","Full_Peak","Region"),transform,max=quantile(Value,1))

regionaldata_Land1.0.1 <- regionaldata_Land1.0a[regionaldata_Land1.0a$ModelN=="MESSAGEix-GLOBIOM_1.0", c(-7)] 
regionaldata_Land1.0.1<- fvarname1(regionaldata_Land1.0.1)

regionaldata_Land1.1a <- regionaldata_Land1.0.1 %>% filter(Year %in% c(2050)) %>% anti_join(Landcoverlist2,by=c("Variable2"))
regionaldata_Land1.1b <- regionaldata_Land1.0.1 %>% filter(Year %in% c(2050)) %>% inner_join(Landcoverlist2,by=c("Variable2"))
regionaldata_Land1.2a <- regionaldata_Land1.0.1 %>% filter(CarbonBudget %in% c(600)) %>% anti_join(Landcoverlist2,by=c("Variable2"))
regionaldata_Land1.2b <- regionaldata_Land1.0.1 %>% filter(CarbonBudget %in% c(600)) %>% inner_join(Landcoverlist2,by=c("Variable2"))

regionaldata_Land2.0 <- regionaldata_Land %>% anti_join(Landcoverlist2,by=c("Variable2"))
regionaldata_Land2.0 <-fprename(regionaldata_Land2.0)
regionaldata_Land2.0 <-modelrename(regionaldata_Land2.0)

g14 <- ggplot() +
  geom_bar(data=subset(regionaldata_Land1.1a,v50>0), aes(x=Region,y=v50,fill=Variable2), position="stack", stat="identity") +
  geom_bar(data=subset(regionaldata_Land1.1a,v50<0), aes(x=Region,y=v50,fill=Variable2), position="stack", stat="identity") +
  geom_errorbar(data=subset(regionaldata_Land1.1b,v50>0), aes(x=Region,ymin=min, ymax=max), colour='gray80', width=.2,position=position_dodge(.9)) + 
  geom_errorbar(data=subset(regionaldata_Land1.1b,v50<0), aes(x=Region,ymin=min, ymax=max), colour='gray80', width=.2,position=position_dodge(.9)) + 
  geom_hline(yintercept=0,linetype="dashed") +
  ylab(expression(paste("million ha"))) + xlab("Region") + MyThemeLine_grid2h +  
  ggtitle(label="Land Use Change in 2050 wrt 2010") +
  theme(legend.position="right", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10)) +
  scale_fill_manual(values=fillpalland)
g14.1 <- g14+facet_wrap(Full_Peak~.)
plot(g14.1)
outname <- paste0("../output/fig/ENGAGE/area/RegionalLUCwrt2010_2050_R5.pdf")
ggsave(g14.1, file=outname, dpi = 600, width=8, height=4,limitsize=FALSE)
outname <- paste0("../output/fig/ENGAGE/pdf/Fig3c.pdf")
ggsave(g14.1, file=outname, dpi = 600, width=8, height=4,limitsize=FALSE)


g15 <- ggplot() +
  geom_bar(data=subset(regionaldata_Land2.0,Value>0), aes(x=ModelN,y=Value,fill=Variable2), position="stack", stat="identity") +
  geom_bar(data=subset(regionaldata_Land2.0,Value<0), aes(x=ModelN,y=Value,fill=Variable2), position="stack", stat="identity") +
  geom_hline(yintercept=0,linetype="dashed") +
  ylab(expression(paste("million ha"))) + xlab("Models") + MyThemeLine_grid2 +
  ggtitle(label="Land Use Change in 2050 with respect to 2010") +
  theme(legend.position="right", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10)) +
  scale_fill_manual(values=fillpalland)
g15.1 <- g15+facet_grid(Full_Peak~Region)
plot(g15.1)
outname <- paste0("../output/fig/ENGAGE/area/RegionalLUCwrt2010_2050_model_R5.png")
ggsave(g15.1, file=outname, dpi = 600, width=12, height=8,limitsize=FALSE)


#----Emis by Region
regionaldata_AFOLU <- filter(worlddata3.0,Variable2 %in% AFOLUlist)
regionaldata_AFOLU <- filter(regionaldata_AFOLU, CarbonBudget %in% c(600), Region!=c("World"), Variable2!=c("Emi_TOT_AFOLU"))
regionaldata_AFOLU$Value <- ifelse(regionaldata_AFOLU$Unit=="Mt CH4/yr", regionaldata_AFOLU$Value*25, regionaldata_AFOLU$Value) 
regionaldata_AFOLU$Value <- ifelse(regionaldata_AFOLU$Unit=="kt N2O/yr", regionaldata_AFOLU$Value*298/1000, regionaldata_AFOLU$Value) 
regionaldata_AFOLU$Unit <- "Mt CO2eq/yr"
regionaldata_AFOLU$Full_Peak <- factor(regionaldata_AFOLU$Full_Peak, levels=c('Peak','Full'))

regionaldata_AFOLU1.0 <- ddply(regionaldata_AFOLU,c("Variable2","Full_Peak","Region","Year"),transform,v50=quantile(Value,0.5))
regionaldata_AFOLU1.0a <- ddply(regionaldata_AFOLU1.0,c("Variable2","Full_Peak","Region","Year"),transform,min=quantile(Value,0))
regionaldata_AFOLU1.0a <- ddply(regionaldata_AFOLU1.0a,c("Variable2","Full_Peak","Region","Year"),transform,max=quantile(Value,1))

regionaldata_AFOLU1.0.1 <- regionaldata_AFOLU1.0a[regionaldata_AFOLU1.0a$ModelN=="MESSAGEix-GLOBIOM_1.0", c(-7)] 
#regionaldata_AFOLU1.0.1 <- regionaldata_AFOLU1.0a[regionaldata_AFOLU1.0a$ModelN=="MESSAGEix-GLOBIOM_1.0",c(1,2,3,4,5,6,7,8,9,10,11,12)] 
regionaldata_AFOLU1.0.1 <- fvarname1(regionaldata_AFOLU1.0.1)

regionaldata_AFOLU1.1a <- regionaldata_AFOLU1.0.1 %>% filter(Year %in% c(2050)) %>% anti_join(AFOLUlist2,by=c("Variable2"))
regionaldata_AFOLU1.1b <- regionaldata_AFOLU1.0.1 %>% filter(Year %in% c(2050)) %>% inner_join(AFOLUlist3,by=c("Variable2"))
regionaldata_AFOLU1.2a <- regionaldata_AFOLU1.0.1 %>% filter(CarbonBudget %in% c(600)) %>% anti_join(AFOLUlist2,by=c("Variable2"))
regionaldata_AFOLU1.2b <- regionaldata_AFOLU1.0.1 %>% filter(CarbonBudget %in% c(600)) %>%inner_join(AFOLUlist3,by=c("Variable2"))
regionaldata_AFOLU1.3 <- filter(regionaldata_AFOLU1.0.1,CarbonBudget %in% c(600),Variable2==c('NET_AFOLU'))

regionaldata_AFOLU <- fvarname1(regionaldata_AFOLU)
regionaldata_AFOLU2.0 <- filter(regionaldata_AFOLU, Year %in% c(2050)) %>% anti_join(AFOLUlist2,by=c("Variable2"))

regionaldata_AFOLU1.2a2<-regionaldata_AFOLU1.2a
regionaldata_AFOLU1.2a$Variable2 <- factor(regionaldata_AFOLU1.2a$Variable2, levels=c('N2O_AFOLU','CH4_AFOLU','CO2_AFOLU','Carbon_Sequestration_BECCS'))
regionaldata_AFOLU1.2a2$Variable2 <- factor(regionaldata_AFOLU1.2a2$Variable2, levels=c('N2O_AFOLU','CH4_AFOLU','Carbon_Sequestration_BECCS','CO2_AFOLU'))

g16 <- ggplot() +
  geom_bar(data=subset(regionaldata_AFOLU1.1a,v50>0), aes(x=Region,y=v50,fill=Variable2), position="stack", stat="identity") +
  geom_bar(data=subset(regionaldata_AFOLU1.1a,v50<0), aes(x=Region,y=v50,fill=Variable2), position="stack", stat="identity") +
  geom_errorbar(data=subset(regionaldata_AFOLU1.1b,v50>0), aes(x=Region,ymin=min, ymax=max), colour='gray80', width=.2,position=position_dodge(.9)) + 
  geom_errorbar(data=subset(regionaldata_AFOLU1.1b,v50<0), aes(x=Region,ymin=min, ymax=max), colour='gray80', width=.2,position=position_dodge(.9)) + 
  geom_hline(yintercept=0,linetype="dashed") +
  ylab(expression(paste("Mt CO2eq/yr"))) + xlab("Region") + MyThemeLine_grid2 +  
  ggtitle(label="Emissions and sequestration in 2050 with respect to 2010") +
  theme(legend.position="right", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10)) +
  scale_fill_manual(values=pastelpal1)
g16.1 <- g16+facet_wrap(Full_Peak~.)
plot(g16.1)
outname <- paste0("../output/fig/ENGAGE/area/RegionalEmiswrt2010_2050_R5.png")
ggsave(g16.1, file=outname, dpi = 600, width=8, height=4,limitsize=FALSE)

regionaldata_AFOLU2.0 <- fprename(regionaldata_AFOLU2.0)
regionaldata_AFOLU2.0 <- modelrename(regionaldata_AFOLU2.0)

g17 <- ggplot() +
  geom_bar(data=subset(regionaldata_AFOLU2.0,Value>0), aes(x=ModelN,y=Value,fill=Variable2), position="stack", stat="identity") +
  geom_bar(data=subset(regionaldata_AFOLU2.0,Value<0), aes(x=ModelN,y=Value,fill=Variable2), position="stack", stat="identity") +
  geom_hline(yintercept=0,linetype="dashed") +
  ylab(expression(paste("Mt CO2eq/yr"))) + xlab("Models") + MyThemeLine_grid2 +
  ggtitle(label="Emissions and sequestration in 2050 with respect to 2010") +
  theme(legend.position="right", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10)) +
  scale_fill_manual(values=pastelpal1)
g17.1 <- g17+facet_grid(Full_Peak~Region)
plot(g17.1)
outname <- paste0("../output/fig/ENGAGE/area/RegionalEmiswrt2010_2050_model_R5.png")
ggsave(g17.1, file=outname, dpi = 600, width=12, height=8,limitsize=FALSE)


# by Time series
g23 <- ggplot() +
  geom_area(data=subset(regionaldata_AFOLU1.2a,v50>0), aes(x=Year,y=v50,fill=Variable2)) +
  geom_area(data=subset(regionaldata_AFOLU1.2a2,v50<0), aes(x=Year,y=v50,fill=Variable2)) +
#  geom_errorbar(data=subset(regionaldata_AFOLU1.0b,v50>0), aes(x=Year,ymin=min/1000, ymax=max/1000), colour='gray80', width=.2,position=position_dodge(.9)) + 
#  geom_errorbar(data=subset(regionaldata_AFOLU1.0b,v50<0), aes(x=Year,ymin=min/1000, ymax=max/1000), colour='gray80', width=.2,position=position_dodge(.9)) + 
#  geom_line(data=regionaldata_AFOLU1.3, aes(x=Year,y=v50), colour="black") +
  geom_hline(yintercept=0,linetype="dashed") +
  ylab("Gt CO2eq/yr") + xlab("Year") + MyThemeLine_grid2 +
  ggtitle(label="GHG Emissions and Sequestration in AFOLU with carbon budget of 600GtCO2") +
  theme(legend.position="right", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10)) +
  scale_fill_manual(values=pastelpal1)
g23.1 <- g23+facet_grid(Region~Full_Peak)
plot(g23.1)
outname <- paste0("../output/fig/ENGAGE/area/RegionalEmission_AFOLU_timeseries.png")
ggsave(g23.1, file=outname, dpi = 600, width=6, height=6,limitsize=FALSE)


#--- Calculate mean value
Landcoverlist_abs_load <- read.table("../define/Landcoverlist_abs.txt", header=F, sep=",") 
Landcoverlist_abs <- as.vector(Landcoverlist_abs_load$V1) 

worlddata_Land_abs <- filter(worlddata3.0,Region %in% c('World'), CarbonBudget %in% c(600), Variable2 %in% Landcoverlist_abs,Year %in% c(2000, 2010, 2020, 2030, 2040, 2050, 2060, 2070, 2080, 2090, 2100))

worlddata_Land_abs <- ddply(worlddata_Land_abs,c("Variable2","Region","Year","Full_Peak"),transform,v50=quantile(Value,0.5))
worlddata_Land_absmean <- worlddata_Land_abs[worlddata_Land_abs$ModelN=="MESSAGEix-GLOBIOM_1.0", c(-1)]

#---


#-- Land difference in FvsP
worlddata_Land3 <- filter(worlddata3.0,Variable2 %in% Landcoverlist)
worlddata_Land3.1 <- filter(worlddata_Land3, CarbonBudget %in% c(600))
worlddata_Land3.1 <- worlddata_Land3.1[c(-2)]
worlddata_Land3.1 <- fvarname1(worlddata_Land3.1)

worlddata_Land3.2 <- worlddata_Land3.1 %>%
  spread(key=Full_Peak, value=Value, fill=0)
worlddata_Land3.2[9] <- worlddata_Land3.2[8]-worlddata_Land3.2[7]
worlddata_Land3.2.1 <- worlddata_Land3.2[c(-7,-8)]

worlddata_Land3.3 <- worlddata_Land3.2.1 %>%
  gather(key=Full_Peak, value=Value,-ModelN,-Region,-Variable2,-Unit,-Year,-CarbonBudget)

#worlddata_Land3.3 <- fvarname1(worlddata_Land3.3)
worlddata_Land3.3.0 <- ddply(worlddata_Land3.3,c("Variable2","Region","Year"),transform,v50=quantile(Value,0.5))
worlddata_Land3.3.0 <- ddply(worlddata_Land3.3.0,c("Variable2","Region","Year"),transform,max=quantile(Value,0.66))
worlddata_Land3.3.0 <- ddply(worlddata_Land3.3.0,c("Variable2","Region","Year"),transform,min=quantile(Value,0.33))
worlddata_Land3.3.1 <- worlddata_Land3.3.0[worlddata_Land3.3.0$ModelN=="MESSAGEix-GLOBIOM_1.0", c(-1)] 
worlddata_Land3.3.2 <- filter(worlddata_Land3.3.1, Year %in% c(2050,2100) & Region!=c("World")) %>% anti_join(Landcoverlist2,by=c("Variable2"))
#worlddata_Land3.3.2 <- filter(worlddata_Land3.3.1, Year %in% c(2050,2100)) %>% anti_join(Landcoverlist2,by=c("Variable2"))
worlddata_Land3.3.3a <- filter(worlddata_Land3.3.1, Region %in% c("World")) %>% anti_join(Landcoverlist2,by=c("Variable2"))
worlddata_Land3.3.3b <- filter(worlddata_Land3.3.1, Region %in% c("World")) %>% inner_join(Landcoverlist2,by=c("Variable2"))

worlddata_Land3.4 <- filter(worlddata_Land3.3, Year %in% c(2050,2100)) %>% anti_join(Landcoverlist2,by=c("Variable2"))

g7 <- ggplot() +
  geom_bar(data=subset(worlddata_Land3.3.2), aes(x=Region,y=v50, fill=Variable2), position="stack", stat="identity") +
  geom_hline(yintercept=0,linetype="dashed") +
  ylab(expression(paste("million ha"))) +
  MyThemeLine_grid2h + 
  ggtitle(label="Land Use Change in Peak relative to Full in 2050") +
  theme(legend.position="right", text=element_text(size=12), axis.title.x = element_blank(),
        legend.title=element_blank()) +
  scale_fill_manual(values=fillpalland)
g7.1 <- g7+facet_grid(.~Year)
plot(g7.1)

outname <- paste0("../output/fig/ENGAGE/area/LUCwrtPvsFbyRegion_R5+world.pdf")
ggsave(g7.1, file=outname, dpi = 600, width=8, height=4,limitsize=FALSE)
outname <- paste0("../output/fig/ENGAGE/pdf/Fib3bRegion.pdf")
ggsave(g7.1, file=outname, dpi = 600, width=8, height=4,limitsize=FALSE)

g8 <- ggplot() +
  geom_bar(data=subset(worlddata_Land3.4), aes(x=ModelN,y=Value, fill=Variable2), position="stack", stat="identity") +
  geom_hline(yintercept=0,linetype="dashed") +
  ylab(expression(paste("million ha"))) +
  MyThemeLine_grid2 + 
  ggtitle(label="Land Use Change in Peak budget relative to Full budget in 2050") +
  theme(legend.position="right", text=element_text(size=12), axis.title.x = element_blank(),
        legend.title=element_blank()) +
  scale_fill_manual(values=fillpalland)
g8.1 <- g8+facet_grid(Year~Region, scales="free_y")
plot(g8.1)
outname <- paste0("../output/fig/ENGAGE/area/LUCwrtPvsFbyRegion_Model_R5.png")
ggsave(g8.1, file=outname, dpi = 600, width=12, height=8,limitsize=FALSE)

# FvsP in Time series


worlddata_Land3.3.3a <- filter(worlddata_Land3.3.3a,Year %in% c(2010,2020,2030,2040,2050,2060,2070,2080,2090,2100))

g21 <- ggplot() +
  geom_bar(data=subset(worlddata_Land3.3.3a,v50>0), aes(x=Year,y=v50,fill=Variable2), position="stack", stat="identity") +
  geom_bar(data=subset(worlddata_Land3.3.3a,v50<0), aes(x=Year,y=v50,fill=Variable2), position="stack", stat="identity") +
#  geom_errorbar(data=subset(worlddata_Land3.3.3b,v50>0), aes(x=Year,ymin=min, ymax=max), colour='gray80', width=.2,position=position_dodge(.9)) + 
#  geom_errorbar(data=subset(worlddata_Land3.3.3b,v50<0), aes(x=Year,ymin=min, ymax=max), colour='gray80', width=.2,position=position_dodge(.9)) + 
  geom_hline(yintercept=0,linetype="dashed") +
  ylab(expression(paste("million ha"))) + xlab("Year") + MyThemeLine_grid2h +
  ggtitle(label="Land Use Change with respect to 2010 with carbon budget of 600GtCO2") +
  theme(legend.position="right", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10)) +
  scale_fill_manual(values=fillpalland)
plot(g21)
outname <- paste0("../output/fig/ENGAGE/area/LUC_timeseriesFvsPv2.pdf")
ggsave(g21, file=outname, dpi = 600, width=6, height=4,limitsize=FALSE)
outname <- paste0("../output/fig/ENGAGE/pdf/Fig3bworld.pdf")
ggsave(g21, file=outname, dpi = 600, width=6, height=4,limitsize=FALSE)

#-- Emissions difference in FvsP

worlddata_AFOLU3 <- filter(worlddata3.0,Variable2 %in% AFOLUlist)
worlddata_AFOLU3$Value <- ifelse(worlddata_AFOLU3$Unit=="Mt CH4/yr", worlddata_AFOLU3$Value*25, worlddata_AFOLU3$Value) 
worlddata_AFOLU3$Value <- ifelse(worlddata_AFOLU3$Unit=="kt N2O/yr", worlddata_AFOLU3$Value*298/1000, worlddata_AFOLU3$Value) 
worlddata_AFOLU3$Unit <- "Mt CO2eq/yr"
worlddata_AFOLU3.1 <- filter(worlddata_AFOLU3, CarbonBudget %in% c(600))
worlddata_AFOLU3.1 <- worlddata_AFOLU3.1[c(-2)]

worlddata_AFOLU3.2 <- worlddata_AFOLU3.1 %>%
  spread(key=Full_Peak, value=Value, fill=0)

worlddata_AFOLU3.2[9] <- worlddata_AFOLU3.2[8]-worlddata_AFOLU3.2[7]
worlddata_AFOLU3.2 <- worlddata_AFOLU3.2[c(-7,-8)]

worlddata_AFOLU3.3 <- worlddata_AFOLU3.2 %>%
  gather(key=Full_Peak, value=Value,-ModelN,-Region,-Variable2,-Unit,-Year,-CarbonBudget)

worlddata_AFOLU3.3 <- fvarname1(worlddata_AFOLU3.3)
worlddata_AFOLU3.3.0 <- ddply(worlddata_AFOLU3.3,c("Variable2","Region","Year"),transform,v50=quantile(Value,0.5))
worlddata_AFOLU3.3.1 <- worlddata_AFOLU3.3.0[worlddata_AFOLU3.3.0$ModelN=="MESSAGEix-GLOBIOM_1.0", c(-1)] 

worlddata_AFOLU3.3.2 <- filter(worlddata_AFOLU3.3.1, Year %in% c(2050,2100) & Region!=c("World")) %>% anti_join(AFOLUlist2,by=c("Variable2"))
worlddata_AFOLU3.3.3a <- filter(worlddata_AFOLU3.3.1, Region %in% c("World")) %>% anti_join(AFOLUlist2,by=c("Variable2"))
worlddata_AFOLU3.3.3b <- filter(worlddata_AFOLU3.3.1, Region %in% c("World")) %>% inner_join(AFOLUlist2,by=c("Variable2"))

worlddata_AFOLU3.4 <- filter(worlddata_AFOLU3.3, Year %in% c(2050,2100)) %>% anti_join(AFOLUlist2,by=c("Variable2"))

g18 <- ggplot() +
  geom_bar(data=subset(worlddata_AFOLU3.3.2), aes(x=Region,y=v50, fill=Variable2), position="stack", stat="identity") +
  geom_hline(yintercept=0,linetype="dashed") +
  ylab(expression(paste("Mt CO2eq/yr"))) +
  MyThemeLine_grid2h + 
  ggtitle(label="AFOLU Emissions in Peak relative to Full in 2050") +
  theme(legend.position="right", text=element_text(size=12), axis.title.x = element_blank(),
        legend.title=element_blank()) +
  scale_fill_manual(values=pastelpal1)
g18.1 <- g18+facet_grid(.~Year)
plot(g18.1)
outname <- paste0("../output/fig/ENGAGE/area/EmiswrtPvsFbyRegion_R5.pdf")
ggsave(g18.1, file=outname, dpi = 600, width=8, height=4,limitsize=FALSE)
outname <- paste0("../output/fig/ENGAGE/pdf/Fig1bRegion.pdf")
ggsave(g18.1, file=outname, dpi = 600, width=8, height=4,limitsize=FALSE)

g19 <- ggplot() +
  geom_bar(data=subset(worlddata_AFOLU3.4), aes(x=ModelN,y=Value, fill=Variable2), position="stack", stat="identity") +
  geom_hline(yintercept=0,linetype="dashed") +
  ylab(expression(paste("Mt CO2eq/yr"))) +
  MyThemeLine_grid2 + 
  ggtitle(label="AFOLU Emissions in Peak budget relative to Full budget in 2050") +
  theme(legend.position="right", text=element_text(size=12), axis.title.x = element_blank(),
        legend.title=element_blank()) +
  scale_fill_manual(values=pastelpal1)
g19.1 <- g19+facet_grid(Region~Year, scales="free_y")
plot(g19.1)
outname <- paste0("../output/fig/ENGAGE/area/EmiswrtPvsFbyRegion_Model_R5.png")
ggsave(g19.1, file=outname, dpi = 600, width=8, height=12,limitsize=FALSE)

worlddata_AFOLU3.3.3a2<-worlddata_AFOLU3.3.3a
#worlddata_AFOLU3.3.3a$Variable2 <- factor(worlddata_AFOLU3.3.3a$Variable2, levels=c('N2O_AFOLU','CH4_AFOLU','CO2_AFOLU','Carbon_Sequestration_BECCS'))
worlddata_AFOLU3.3.3a2$Variable2 <- factor(worlddata_AFOLU3.3.3a2$Variable2, levels=c('N2O_AFOLU','CH4_AFOLU','CO2_AFOLU','Carbon_Sequestration_BECCS'))


# FvsP in Time series
g22 <- ggplot() +
  geom_area(data=subset(worlddata_AFOLU3.3.3a,v50>0), aes(x=Year,y=v50,fill=Variable2)) +
  geom_area(data=subset(worlddata_AFOLU3.3.3a,v50<0), aes(x=Year,y=v50,fill=Variable2)) +
  #  geom_errorbar(data=subset(worlddata_AFOLU3.3.3b,v50>0), aes(x=Year,ymin=min, ymax=max), colour='gray80', width=.2,position=position_dodge(.9)) + 
  #  geom_errorbar(data=subset(worlddata_AFOLU3.3.3b,v50<0), aes(x=Year,ymin=min, ymax=max), colour='gray80', width=.2,position=position_dodge(.9)) + 
  geom_hline(yintercept=0,linetype="dashed") +
  ylab("Gt CO2eq/yr") + xlab("Year") + MyThemeLine_grid2h +
  ggtitle(label="GHG Emissions and Sequestration in AFOLU with carbon budget of 600GtCO2") +
  theme(legend.position="right", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10)) +
  scale_fill_manual(values=pastelpal1)
plot(g22)
outname <- paste0("../output/fig/ENGAGE/area/Emission_AFOLU_timeseriesFvsP.pdf")
ggsave(g22, file=outname, dpi = 600, width=6, height=4,limitsize=FALSE)
outname <- paste0("../output/fig/ENGAGE/pdf/Fig1bworld.pdf")
ggsave(g22, file=outname, dpi = 600, width=6, height=4,limitsize=FALSE)


#--Food analysis

Analizedvarlist_load <- read.table("../define/Analizedvarlist.txt", header=F, sep=",") 
Analizedvarlist <- as.vector(Analizedvarlist_load$V1) 

worlddata_Food <- filter(worlddata3.0,Region %in% c('World'),Variable2 %in% Analizedvarlist,Year %in% c(2000, 2010, 2020, 2030, 2040, 2050, 2060, 2070, 2080, 2090, 2100))
#worlddata_Food <- fvarname1(worlddata_Food)

worlddata_Food0.1 <- filter(worlddata_Food,Year %in% c(2050))
worlddata_Food0.2 <- filter(worlddata_Food,CarbonBudget %in% c(600))

worlddata_Food1.0 <- ddply(worlddata_Food,c("Variable2","Full_Peak","CarbonBudget","Year"),transform,v50=quantile(Value,0.5))
worlddata_Food1.0.1 <- worlddata_Food1.0[worlddata_Food1.0$ModelN=="AIM_CGE 2_2", c(-7)] 
worlddata_Food1.1 <- filter(worlddata_Food1.0.1,Year %in% c("2050","2100"),CarbonBudget %in% c(600))
worlddata_Food1.2 <- filter(worlddata_Food1.0.1,CarbonBudget %in% c(600))

worlddata_Food1.1 <- worlddata_Food1.1[c(-2)]
worlddata_Food1.1 <- worlddata_Food1.1 %>%
  spread(key=Full_Peak, value=v50, fill=0)

worlddata_Food1.1[9] <- worlddata_Food1.1[8]-worlddata_Food1.1[7]
worlddata_Food1.1 <- worlddata_Food1.1[c(-8)]
names(worlddata_Food1.1) <- c("ModelN","Region","Variable2","Unit","Year","CarbonBudget","Full","Changes in Peak from Full")	#Rename the column
worlddata_Food1.1 <- worlddata_Food1.1 %>%
  gather(key=Full_Peak, value=Value,-ModelN,-Region,-Variable2,-Unit,-Year,-CarbonBudget)

worlddata_Food1.1$Full_Peak <- factor(worlddata_Food1.1$Full_Peak,levels=c('Full','Changes in Peak from Full'))
worlddata_Food1.1 <- fvarname1(worlddata_Food1.1)
worlddata_Food1.1$Year <- as.character(worlddata_Food1.1$Year)

linepalette_FP=c('Full'='red','Peak'='blue','Changes in Peak from Full'='blue')

g20 <- ggplot() +
  geom_bar(data=subset(worlddata_Food1.1), aes(x=Year,y=Value,fill=Full_Peak),position="stack", stat="identity") +
  geom_hline(yintercept=0,linetype="dashed") +
  ylab("") + xlab("Year") + MyThemeLine_grid2 +  
  ggtitle(label="Effects on Food systems") +
  theme(legend.position="bottom", text=element_text(size=12),legend.title=element_blank(),legend.text = element_text(size = 10),
        axis.text.x=element_text(angle=0, vjust=0.5, hjust=0.5)) +
  scale_fill_manual(values=linepalette_FP)
g20.1 <- g20+facet_grid(Full_Peak~Variable2, scales="free_y")
plot(g20.1)
outname <- paste0("../output/fig/ENGAGE/bar/FoodSec_20502100.png")
ggsave(g20.1, file=outname, dpi = 600, width=8, height=6,limitsize=FALSE)


#----




worlddata1.0 <-worlddata %>%
  rename("2000"=X2000,"2005"=X2005,"2010"=X2010,"2015"=X2015,"2020"=X2020,"2025"=X2025,"2030"=X2030,"2035"=X2035,"2040"=X2040,"2045"=X2045,"2050"=X2050,"2055"=X2055,"2060"=X2060,"2065"=X2065,"2070"=X2070,"2075"=X2075,"2080"=X2080,"2085"=X2085,"2090"=X2090,"2095"=X2095,"2100"=X2100) %>% 
  gather(key=Year,value=Value,-Model,-Scenario,-Region,-Variable,-Unit) %>% inner_join(variablemap,by="Variable")  %>%
  select(-Variable) %>% inner_join(modelemap,by="Model") %>% select(-Model)  %>% select(append(colhead2,c("Year","Value")))
worlddata2.2 <- filter(worlddata2.1,Year %in% c(2010,2030,2050,2100))
